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ABSTRACT 

A  numerical  code,  USPOTF2,  has  been  formulated  to  solve  for  the  potential  flow 
for  two  airfoils  executing  unsteady  motions  in  an  inviscid  incompressible  flow  medium. 
This  code  is  an  extension  of  an  existing  code  U2DIIF,  which  does  the  same  calculations 
for  the  single  airfoil  case.  The  technique  uses  the  well  known  Panel  Methods  for  steady 
flow  and  extends  it  to  unsteady  flow  by  introducing  a  wake  model  which  creates  a  non- 
linear problem  due  to  the  continuous  shedding  of  vortices  into  the  trailing  wake.  The 
presence  of  the  second  airfoil  introduces  a  set  of  non-linear  coupled  equations  for  the 
Kutta  condition.  Numerous  case-runs  are  presented  to  illustrate  the  capability  of  the 
code.  The  case  of  the  step  change  in  angle  of  attack  is  compared  with  Giesing's  work. 
All  other  case-runs  are  illustrated  together  with  the  results  for  the  single  airfoil  case. 
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I.     INTRODUCTION 

A.  GENERAL 

In  this  paper  ,  a  numerical  method  is  formulated  to  solve  for  the  flow  about  two 
two-dimensional  airfoils  which  are  arbitrary  located  and  are  performing  an  arbitrary  time 
dependent  motion  in  an  inviscid  incompressible  fluid.  The  original  work  by  Teng  [Ref. 
1]  is  for  a  single  arbitrarily  defined  airfoil  in  the  same  potential  flow  condition.  The  ex- 
tension of  Teng's  code  to  two  bodies  is  considered  here.  Where  possible,  the  same  con- 
vention and  notation  are  adopted.  As  for  the  single  airfoil  case,  all  velocities  are 
non-dimensionalised  with  respect  to  the  free  stream  and  all  lengths  with  respect  to  the 
chord  length. 

A  new  subroutine  (SUBROUTINE  NEWPOS)  is  added  to  transform  either  of  the 
two  local  coordinate  systems  to  the  global  coordinate  system.  A  modification  is  entered 
into  the  treatment  of  the  Kutta  condition  to  make  it  consistent  with  the  unsteady  flow 
Kutta  condition  treatment.  Subroutine  COF1SH  is  deleted  and  its  role  is  taken  over  by 
Subroutine  COEF  which  normally  performs  the  formulation  for  the  flow  tangency  con- 
dition for  the  unsteady  case.  A  more  accurate  method  is  also  introduced  to  obtain  the 
velocity  potential.  The  reader  is  referred  to  the  work  of  Krainer  [Ref.  2j  for  further 
improvement  of  the  original  code. 

This  documentation  is  set  up  as  for  the  original  documentation  in  order  that  it  will 
be  easier  to  follow  and  cross-referenced.  While  it  is  the  intent  of  the  author  to  keep  this 
thesis  as  complete  as  possible,  the  reader  is  well  advised  to  review  Reference  1  for  the 
work  involved  in  the  single  airfoil  case;  no  special  effort  will  be  made  to  reproduce  it 
here. 

B.  BASIC  THEORY  AND  APPROACH 

1.     Steady  Flow  Problem 

The  treatment  for  the  two  airfoils  case  in  steady  incompressible  flow  follows 
closely  to  the  single  airfoil  case.  The  governing  equation,  as  for  the  single  airfoil  case  , 
follows  from  the  Conservation  of  Mass  and  the   Condition  of  Irrotational  Flow. 

The  continuity  equation  of  an  incompressible  fluid  (div  V  =  0)  and  the  condi- 
tion of  irrotational  flow  (curl  V  =  0)  leads  to  the  well  known  Laplace  Equation. 

div(grad<?  j  =  0  (1.1) 


with  </>  denoting  the  disturbance  potential  for  the  velocity.  This  is  seen  to  be  the  classic 
Neumann  problem  of  potential  theory  with  the  usual  problem  of  defining  the  boundary 
conditions.  The  boundary  conditions  for  the  disturbance  potential  <£  are  that  its  gradi- 
ent normal  to  the  surface  be  equal  to  the  normal  velocity  of  the  surface  and  that  its 
gradient  vanish  at  infinity;  that  is 

V0./7  =  Vr.n  onS  (1.2) 

lim  V<t>{P)  -  0  (1.3) 

P->  CO 

where  V,  is  the  resultant  velocity  of  a  point  on  the  body  as  seen  from  the  inertia  frame 
of  reference,  n  is  the  outward  unit  vector  normal  to  the  body,  S  is  the  body  surface  and 
P  represents  a  general  point.    Equation  1.2  holds  for  both  airfoils  i.e.  on  both  surfaces. 

The  pressure  coefficient  is  obtained  through  the  Bernoulli's  equation  which  de- 
rives from  the  Momentum  equation. 

The  approach  adopted  is  associated  with  Hess  and  Smith  [Ref.  3]  who  devised 
the  popularly  known  PAS  EL  method  in  the  early  sixties.  In  words,  the  boundary  or 
airfoil  surfaces  S1  and  S2  about  which  the  flow  is  to  be  computed  is  approximated  by  a 
large  number  of  surface  elements  whose  characteristic  dimensions  are  small  compared 
to  those  of  the  body.  Over  each  surface  element,  a  uniform  source  distribution  and  a 
uniform  vorticity  distribution  is  placed.  The  source  strength  (q.)  varies  from  element  to 
element,  while  the  vortex  strength  {■■;,)  is  the  same  for  all  elements  in  the  same  airfoil  but 
is  different  across  the  airfoil.  The  singularity  strengths  are  determined  from  the  flow 
tangency  condition  on  both  body  surfaces  and  the  two  Kutta  conditions  at  both  trailing 
edges.  With  the  determination  of  the  singularity  strengths,  the  relevant  aerodynamic 
data  can  then  be  subsequently  computed. 
2.     Unsteady  Flow  Problem 

The  unsteady  problem  is  similar  to  the  steady  flow  problem  in  that  they  both 
have  the  same  governing  equation  viz.  the  Laplace  equation,  and  that  for  both  problems, 
the  pressure  and  velocity  are  decoupled  so  that  the  velocity  and  pressure  calculation  can 
be  computed  separately  and  consecutively. 

This  problem  differs  from  the  steady  flow  in  that  another  model  is  required  to 
simulate  the  continuous  shedding  of  vorticity  into  the  trailing  wake.  The  existence  of  a 
vortex  sheet  behind  the  airfoil  can  be  explained  by  the  Helmholtz  theorem  which  is 
basically  a  statement  of  the  Conservation  of  Vorticity.   This  requires  that  any  change  in 


the  circulation  around  the  airfoil  must  be  matched  by  an  equal  and  opposite  vortex 
somewhere  in  the  flowfield.  The  presence  of  the  countervortices  provides  the  flow  with 
a  kind  of  a  memory  in  that  the  flow  at  a  particular  time  is  affected  by  the  bound  circu- 
lation of  the  past.  It  is  this  non-linearity  that  distinguishes  the  numerical  technique  from 
the  simple  steady  flow  problem  of  solving  N  linear  equations  in  N  unknowns. 

The  solution  technique  requires  an  iterative  type  solution.  The  present  ap- 
proach follows  closely  the  original  panel  method  of  Hess  and  Smith  as  described  in  the 
steady  flow  development,  while  with  regard  to  the  modelling  of  the  wake,  it  adopts  the 
procedure  advocated  by  Basu  and  Hancock  [Ref.  4J.  A  uniform  source  distribution  (q,)k 
and  a  uniform  vorticity  distribution  [>•(/)]*  as  for  steady  flow  is  placed  on  each  panel  at 
time  t,  where  j  denotes  the  panel  number  and  /  the  airfoil  number.  The  wake  consists 
of  a  single  vorticity  panel  attached  as  an  additional  element  on  each  airfoil  through 
which  the  vorticities  are  shed  into  the  respective  wake  and  a  series  of  point  vortices 
which  are  being  convected  downstream  with  the  fluid.  A  uniform  vorticity  distribution 
of  strength  (/J/)),  is  placed  on  the  wake  panel  of  each  airfoil.  This  panel  is  further 
characterised  by  its  length  A{f)k  and  its  inclination  (  0(„-d,  )*  with  respect  to  the  respective 
local  frame  of  reference.  After  each  time  step,  the  vorticity  of  the  wake  panel  is  con- 
centrated into  a  single  point  vortex  and  convected  downstream.  Simultaneously  a  new 
wake  panel  is  formed.  The  downstream  wake  of  point  vortices  is  thus  formed  by  the 
shed  vorticity  of  previous  time  steps. 

C.     SCOPE 

Chapter  II  extends  the  original  code  to  handle  the  steady  flow  problem  for  two 
airfoils  set  at  different  relative  distances  and  angles  of  attack. 

Chapter  III  deals  with  the  unsteady  problem  for  the  two  airfoils  system.  It  intro- 
duces a  new  subroutine  and  a  more  accurate  method  of  calculating  the  velocity  poten- 
tial. The  Kutta  condition  for  the  two  airfoils  system  is  specially  treated  as  its  unique 
problem  of  a  non-linear  coupled  system  is  not  seen  in  the  single  airfoil  case. 

Chapter  IV  describes  the  computer  program,  its  essential  capabilities  and  limita- 
tions, its  associated  subroutines,  its  input  requirements  and  its  associated  output  print- 
out. 

Chapter  V  presents  the  results  of  some  case-runs.  Of  interest  is  a  comparison  case 
of  a  step  change  in  angle  of  attack  (AOA)  with  Giesing  [Ref.  5]  for  the  same  airfoil 
undergoing  an  impulsive  start  at  the  same  AOA.  Case-runs  will  also  be  run  with  both 
airfoils  at  large  distances  apart  to  compare  with  the  single  airfoil  case.    In  addition,  ex- 


ample  cases  for  which  no  comparisons  exist  are  given,  to  exhibit  the  capability  of  the 
method. 

Finally,  Chapter  VI  concludes  with  future  development  efforts  and  the  application 
potential  of  this  numerical  method. 


II.     STEADY  FLOW  PROBLEM  FORMULATION  FOR  TWO  AIRFOILS 

A.  GENERAL 

The  modification  work  for  the  steady  flow  is  straightforward.  The  revised  program 
allows  for  two  arbitrarily  defined  airfoils  placed  at  an  arbitrary  distance  set  at  different 
angles  of  attack.  For  reason  of  simplicity,  the  number  of  panels  and  nodes  and  the  pivot 
location  are  set  to  be  the  same  for  both  airfoils. 

B.  FRAMES  OF  REFERENCE 

Three  frames  of  reference  are  involved  in  the  two  two-dimensional  airfoils'  case  in 
steady  flow.   These  are  three  inertia  frames  of  references  as  indicated  in  Figure  1. 

The  first  inertia  frame  of  reference  (also  known  as  global  frame  of  reference)  is  set 
at  the  pivot  position  of  the  first  airfoil  with  the  X-axis  pointing  in  the  direction  of  the 
free-stream  velocity.  The  two  other  inertia  frames  of  references,  henceforth,  will  be 
known  as  two  frozen  local  frames  of  referencel  (xKv,  and  x.  v2)  are  fixed  respectively  to 
each  airfoil  with  the  x-a\is  coinciding  with  the  chord  line  onginating  from  the  respective 
leading  edge.  The  two  local  frames  of  reference  are  set  apart  by  XShift  and  YShift  on 
the  global  frame  of  reference.  In  steady  flow,  the  fluid  velocity  and  pressure  depend  only 
on  the  spatial  coordinates  (X,Y)  and  not  on  time. 

The  two  airfoils  are  defined  in  the  local  coordinate  system  as  input  data  for  sim- 
plicity and  are  then  transformed  to  the  global  coordinate  system  through  a  knowledge 
of  the  relative  positions  of  the  2  airfoils'  pivots  positions  and  the  respective  local  angles 
of  attack. 

C.  STEADY  FLOW  PANEL  METHODS 
I.     Definition  of  nodes  and  panels 

Each  airfoil  surface  is  divided  into  n  straight  line  segments  called  panels  by 
(n+1)  arbitrarily  chosen  points  called  nodes.  The  numbering  sequence  begins  with 
panel  1  on  the  lower  surface  at  the  first  airfoil  trailing  edge  and  proceeds  clockwise 
around  the  airfoil  contour  so  that  the  last  panel  on  airfoil  1  ends  on  the  upper  surface 
at  the  trailing  edge.  This  numbering  sequence  then  proceeds  in  a  similar  fashion  for  the 
second  airfoil  (  See  Figure  2).    As  with  the  single  airfoil  case,  the  numbering  sequence 


1  For  the  steady  case,  the  notation  frozen'  will  be  dropped 


Figure  1.      Frames  of  Reference  for  Steady  Flow. 


dictates  that  the  airfoil  body  always  lies  on  the  right  hand  side  of  the  i"'  panel  as  one 
proceeds  from  the  i*  node  to  the  (i+  1)"'  Also  the  1*  and  {n  +  l)r/'  nodes  coincide  at  the 
first  airfoil  trailing  edge  and  the  {n  +  2)!h  and  {In  +  2)"1  coincide  at  the  second  airfoil 
trailing  edge.  This  numbering  system  facilitates  the  definitions  of  the  unit  normal  vector 
n,  and  the  unit  tangential  vector  t,  for  all  panels  with  n,  being  directed  outward  from  the 
body  into  the  flow  and  /,  directed  from  the  i'h  node  to  the  (/+  l)'h  node.  The  numbering 
of  the  panel  system  is  somewhat  complicated  by  the  fact  that  a  continuous  panel  num- 
bering sequence  across  the  airfoil  is  desired.  This  procedure  leads  to  the  peculiar  panel 
numbering  behaviour  seen  in  Figure  2  where  for  the  first  airfoil,  the  /'"'  panel  lies  between 
/'"■  and  (i  +  l)r*  nodes  while  for  the  second  airfoil,  the  i'h  panel  lies  between 
(/  +  1)*  and  (/  +  2)"'  nodes-. 

2.     Distribution  of  Singularities 

Over  each  surface  element  of  the  two  airfoils,  a  uniform  source  distribition  and 
a  uniform  vorticity  distribution  is  placed.  The  source  strength  q,  varies  from  element  to 
element  for  each  airfoil  while  the  vorticity  strength  y,  remains  the  same  for  all  elements 
in  the  same  airfoil  but  is  different  across  the  airfoil.  This  choice  of  singularities  follows 
closely  the  original  panel  method  of  Hess  and  Smith.  It  automatically  satisfies  the 
Laplace  Equation  (which  is  the  governing  equation  for  the  inviscid  incompressible  flow) 
and  the  boundary  condition  at  the  far  field  (oo).  In  addition,  as  the  Laplace  Equation 
is  a  linear  homogeneous  second  order  partial  differential  equation,  an  overall  compli- 
cated flow  field  can  be  built  up  by  the  combination  of  simple  flows  with  the  condition 
that  the  appropriate  boundary  condition  on  the  airfoil  be  satisfied  accurately. 

For  our  case,  the  overall  flow  field  (represented  by  the  velocity  potential  O)  can 
be  built  up  by  three  simple  flows.  Writing  this  in  terms  of  the  respective  local  frame  of 
reference. 

<to*  .y)  =  4>J.x  ,y)  +  <t>s(*  .jO  +  4>A*  .jO  (2-j) 

where  4>x{x  ,  y)  is  the  potential  of  the  onset  flow, 

*«,(*  •>)  =  rjfiO*  «(0  +ySm  «(/)]  (2-2) 


2  For  the  code  8,  has  been  defined  out  of  the  above  convention  with  d„  i=  l,2...n  as  per  panel 
numbers  for  the  first  airfoil,  but  with  0„_,  reserved  for  the  wake  element  of  the  first  airfoil  and  6„ 
i  =  n+  2.n  +  3...2n+  1  for  the  second  airfoil  with  0^_2  reserved  for  the  wake  element  of  the  second 
airfoil  for  unsteadv  flow. 


J+1 


1 


Panel  j 

Source  distribution  =  q. 

Vorticity  distribution  =  7. 


Note:  1. Nodal  points  defined  by  i  =  1,2 

and  i  =  n+2 


n+1  for  first  airfoil 
2n+2  for  second  airfoil. 


2. Panel  number  defined  by  j  =  1,2  ...  n  for  first  airfoil 

and  j  =  n+1  ...  2n  for  second  airfoil. 
Figure  2.   Panel  Methods  Representation  for  Steady  Flow. 


4>s  is  the  velocity  potential  of  a  source  of  distribution  q(s)  per  unit  length. 


J 


—  Inrds  (2.3) 


<f>,  is  the  velocity  potential  of  a  vorticity  distribution  y  (s)  per  unit  length 


-• 


4>,  =  - 


—  e  ds  (2.4) 


At  this  point,  the  disturbance  potential,  </>,  is  introduced,  which  is  defined  to  be  the  sum 
of  the  potential  due  to  the  source  and  vorticity  distribution. 

Q  =  (t>s  +  4>v  (2.5) 

Equation  (2.1)  can  then  be  read  as 

4>  =  0x  +  4>  (2.6) 

The  convenience  of  defining  the  above  allows  for  the  total  velocity  vector  to  be  viewed 
as  two  components  viz.  the  onset  and  the  induced  velocity  due  to  the  disturbance  po- 
tential.  The  total  velocity  is  thus  : 

=  V6oc  +  V^>  (2.7) 

It  is  in  the  introduction  of  the  disturbance  potential  that  leads  us  to  the  concept  of  in- 
fluence coefficient  which  will  be  elaborated  in  a  later  section. 

The  pressure  coefficient  can  be  obtained  from  Bernoulli's  Equation  which  is 
derived  from  the  Conservation  of  Momentum. 


JHzIr 

1     i/ 


c°=1r-7r  =  x-\~vf}2  (2-8) 


3.     Boundary  Condition 

As  in  the  single  airfoil  case,  the  boundary'  conditions  to  be  satisfied  include  the 
flow  tangency  conditions  and  the  Kutta  condition.  The  flow  tangency  conditions  are 
satisfied  at  the  exterior  mid-points  (control  points).  The  normal  velocity  is  taken  with 
respect  to  the  respective  local  frame  of  reference  (for  consistency  with  unsteady  flow 
notation  to  be  introduced  later): 

(Vn);  =  0,      /=  1,2  ...«,«+  1  ...2«  (2.9) 

where  ( Vn),  is  the  normal  component  of  the  total  velocity^. 

The  Kutta  condition  postulates  that  the  pressure  on  the  upper  and  lower  sur- 
face at  the  trailing  edge  of  each  airfoil  be  equal.  For  steady  potential  flow,  equal  pres- 
sure implies  equal  tangential  velocity  in  the  downstream  direction  at  the  first  and  last 
panel  of  each  airfoil  viz  the  Bernoulli  equation.  With  our  definition  of  the  tangential 
vector  we  then  have 

{V\    =-(V')n  1  si  airfoil 

{Vl)n^  =  -{V\n       2nd  airfoil  (2.10) 

As  with  the  single  airfoil,  equations  2.9,  and  2.10  lead  to  a  linear  system  of 
(2n+  2)  simultaneous  equations.   With  Vn  and  V'  expressed  explicitly  in  terms  of  q,  (j  = 
1,  2, ...  n,  n+  1...  2n)  and  y,(l=  1,  2)  we  have  2n  +  2  unknowns  in  (2n  +  2)  system  of  linear 
simultaneous  equations  which  can  be  easily  solved. 

D.     INFLUENCE  COEFFICIENT 

1.     Concept  of  Influence  Coefficient 

As  introduced  earlier,  this  important  concept  of  influence  coefficient  results 
from  the  presence  of  the  disturbance  potential  which  follows  from  the  presence  of  the 
singularities.  Formally,  an  influence  coefficient  is  defined  as  the  velocity  induced  at  a 
field  point  by  a  unit  strength  singularity  (be  it  a  point  singularity  or  a  distributed 
singularity)  placed  anywhere  within  the  flow  field.  Recall  that  equations  2.9  and  2.10 
require  the  computation  of  the  normal  and  tangential  velocity  components  at  all  the  el- 
ements control  points.    The  normal  components  of  velocities  are  essential  in  satisfying 


3  This  follows  Giesing's  notation  where  the  velocity  with  respect  to  the  airfoil  frame  is  denoted 
as  total  velocity.  While  there  is  no  misunderstanding  in  the  steady  flow,  the  notation  actually  refers 
to  the  relative  velocity  with  respect  to  the  airfoil  moving  frame  of  reference  for  the  unsteady  flow. 
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the  flow  tangency  conditions,  while  the  tangential  components  of  velocities  are  necessary 
for  satisfying  the  Kutta  condition  as  well  as  computing  the  pressure  distribution. 

2.  Notation  for  Influence  Coefficient 

The  notation  used  in  this  documentation  will  be  as  for  the  notation  for  the 
single  airfoil  case  with  a  slight  modification  to  take  into  account  the  effects  of  the  the 
second  airfoil.  As  the  influence  coefficients  are  related  to  the  geometry  of  the  airfoil  and 
their  relative  positions,  it  must,  of  necessity  be  computed  with  respect  to  the  global 
frame  of  reference  for  the  two  airfoils'  case. 

For  steady  flow,  tne  following  influence  coefficients  are  defined. 

•  A".  :    normal  velocity  component  induced  at  the  /"'  control  point  by  unit  strength 
source  distribution  on  the/'-  panel. 

•  A;  :  tangential  velocity  component  induced  at  the  i'h  control  point  by  unit  strength 
source  distribution  on  the/''  panel. 

•  Z?".  :    normal  velocity  component  induced  at  the  i'1'  control  point  by  unit  strength 
vorticity  distribution  on  the/*  panel. 

•  B"  :  tangential  velocity  component  induced  at  the  /'*  control  point  by  unit  strength 
vorticity  distribution  on  the/''  panel. 

where  1  and  j  denotes  panel  numbers  and  are  defined  as  : 
r  =  I,  2, ...  n,  n+1  ...  2n 

/'■  =   1.  2.  ...  n.  n+1  ...  2n 

3.  Computation  of  Influence  Coefficient 

A  single  source  located  on  the/'1  panel  in  the  local  frame  of  reference  (see  Figure 
3)  induces  a  total  velocity  of 

\-=^-  (2.11) 

where  the  component  perpendicular  to  the  inducing  panel  is 

and  the  component  parallel  to  the  inducing  panel  is 

2nr  '  In      h-  +  52 

For  a  distributed  source  panel  on  the/''  panel  in  the  local  frame  of  reference,  we  have 
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n; 


-4 


(Xj  •  Yj) 


jth  panel 


2  J 


(Xj+1   .  yj+O 


Figure  3.      Influence  Coefficient  for  Source  Panel 
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Vj    =    \\{s)ds  =^  FTAX  (2.14) 

where 

_,  (-xmi  ~  Xj){ymi  ->}+,)  -  (ym,  - yj){xmt  -  xJ+l) 

FT  AS    =  tan : : - 

{xmi  -  Xj){xml  -  xy+1)  +  (jv?j(-  -  yfjiym-,  -yJ+l) 

xm,  =  l/2(xf-hxf+1) 

ym,  =  I/2fo+Jfo.i)  (2.15) 


Kis)^  =    -T7-  FLOG  (2.16) 


1    -  f    ^ 


->_ 


where 

FLOG    =    In ! ^  , ^-  (2.17) 

(jwnsi  —  jcy)  +(yml-yj) 

Transforming  the  above  to  global  frame  of  reference,  we  have 

A%  =   VjCoM-Bjl  -  rJSin(0,-0;)  (2.18) 

4  =   IpmiOt-ej)  +  VfioiQt-Oj)  (2.19) 

^  =    -  ( I  TjCosidi  -  6j)  -  I  {Shift  -  6j))  (2.20) 

B[j  =   PjSinft  -  6j)  -  KjCosft  -  6j)  (2.21) 

For  the  Steady  Flow  code,  we  define 

AAX(IJ)  =  A-j  =  By  (2-22) 

BB.X(IJ)  =    -A[j  =  By  (2-23) 


13 


E.     NUMERICAL  SOLUTION  SCHEME 

1.     Rewriting  the  Boundary  Condition 

Though  it  is  not  critical  for  the  steady  flow  problem  to  adopt  any  particular 
frame  of  reference,  we  will  adopt  the  local  frame  of  reference  in  satisfying  the  boundary 
condition  for  consistency  with  the  treatment  used  in  the  unsteady  case. 

Even  though  the  influence  coefficients  are  defined  in  terms  of  the  global  frame 
of  reference,  the  computation  of  the  normal  and  tangential  velocities  of  each  panel  is 
independent  of  the  coordinate  system  used.  Thus  the  normal  and  tangential  velocities 
obtained  with  the  global  coordinates  system  would  be  the  same  as  those  obtained  for  the 
local  coordinate  system.  Using  the  local  coordinate  system,  the  flow  tangency  condition 
of  equation  2.9  is  defined  as  follows  : 

n  In  n  In 

{V\  =    jMjfa  +   Yj  AM  +  V(1)X^  +  K2)  X  ^  +  V°°  Sln[a(/)  "  °i]   =   0        (2"24) 

where 

/'=  1,2  ...  «      ;  /  =  1, 

=  n+  1  ...  In  ;     =2 

The  Kutta  condition  of  eqn  2.10,  in  terms  of  influence  coefficients,  becomes  for  airfoil 
1: 

n  In  n  In 

£j\flt  +  2]  A\A  +  ywYjBv +  y(2)  X  B'j  +  v°° Cos[a(1)  ~  e>]  = 

y=1  y=n+1  j=\  j=n+\ 

n  In  n  In 

-  \Yja'm  +  Yj  A'njqj + y(i)ZX + y(2)  Yj  B'nj + v°° cos[a(i)  ~  dni\     (2,25) 

W=l  j=n+\  j=\  j=n+\  J 

and  for  airfoil  2: 

n  2n  n  2n 

Yja"+>m+  X^+i^+y(i)S^+,^+y(2)Z^+i^+ro°cos[a(2)~^+i]  = 

j=\  j=n+l  j=\  j=n+] 
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_f 


2n 


j  2/j./g +  X!  ^^ + ^z/2'^ + >,(2)  X  ^ + r~  c°s[a(2)  -  ^ j  (2-26) 

2.     Solving  for  the  Strengths  of  Source  and  Vorticity  Distribution 

Equations  2.24,  2.25  and  2.26  can  be  written  as  a  set  of  2n  +  2  linear  simultane- 
ous equations  with  2n+2  unknowns  {qj,j=  1,2  ...  2A?andyy,/  =  1,2)  and  solved  as  for 
the  single  airfoil  case.  However,  to  make  the  routine  consistent  with  the  unsteady  flow 
case,  the  following  method  is  adopted. 

The  flow  tangency  condition  can  be  rewritten  explicitly  for  the  q,  in  terms  of 
y(l),  y(2)  and  the  free  stream  constant  term;  that  is 


au     a, 


'2.1 


1.2 
a2,2 


a 


3.2 


42.3 


*3,3 


l2n,\     u2n,2     u2n,3 
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- 

_ 
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(2.27) 

Gauss  Elimination  is  then  used  to  solve  for  the  qs  in  terms  of  y(l),  y(2)  and  the 
constant  term.  This  eives 


qj  =  b\jy{\)  +  b2jV(2)  +  blj  j  -    1,2 ...  2« 


(2.28) 


Equation  2.28  is  then  substituted  into  the  Kutta  condition  at  the  two  trailing  edges  to 
form  two  linear  simultaneous  equation  with  two  unknowns  y(l)andy(2)  .  Gauss 
Eimination  is  again  used  to  solve  for  the  Vorticity  distribution  y(l)andy(2)  and  these 
are  then  back  substituted  to  solve  for  the  q.. 

3.     Computation  of  Velocity  and  Pressure  Distribution 

Once  the  q.  (J  =  1.2  ...  2/?)  and  y,  (I  =  1,2)  are  solved,  the  velocities  at  all  the  panel 
control  points  can  be  easily  obtained.    The  normal  velocity  is  given  by 
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n  In  n  2n 

{V\  =  Y/k  +  Yj  A"}q> +  y(1)X^  +  y(2)  X  ^  +  *  ~ Sin[  a(/) ~  6i  ]  (2'29) 

while  the  tangential  velocity  is  given  by 

n  In  n  In 

[y\  =  YAk+  YA'M+y{1)YjB'u*y{2)YjB'J+  v°oC°sl«(r)-ei)         (2-30) 

y=l  j=n+\  j=\  j=n+\ 

with 

/=  1,2...  n     ;  /  =  1, 
=  «  +  1  ...  In  ;     =2 

The  normal  velocity  will  obviously  satisfy  equation  2.9  (used  to  check  the  code)  showing 
that  the  flow  tangency  condition  is  satisfied  while  the  the  total  velocity  will  be  given  by 
the  tangential  velocity. 

Vtotai   =  {Vl)t,    i=\2...2n  (2.31) 

where 


{V)t  =  -  YjBm  "  Yj  Bvqj  +  y(1)Z^  + 1/(2)  X  AS +  K~Cos[  a(/)  ~ e' ]      (2-32) 

Substituting  equation  2.31  into  2.8  for  C„  with  (V1),  defined  as  in  equation  2.32,  the 
pressure  coefficient  at  the  i'h  control  point  is 

{Cp)t  =  \-{V1)],    /-1.2...2/I  (2.33) 

4.     Computation  of  Forces  and  Moments 

The  two-dimensional  aerodynamic  lift  (Q,  drag  (CJ  and  pitching  moment  (CJ 
are  calculated  with  respect  to  the  global  frame  of  reference.  The  moment  coefficients 
are  computed  with  respect  to  the  respective  leading  edges.   For  the  code,  we  have 

n 

c/0  =  YscpMxm  -  Xi)  (2-34) 

;=i 
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QO  =    -X(C^O;+1-r()  (2.35) 

(=1 

n 

Cm(f)  =  JjCJA  (Xi+l  -  A-,.)Am,  +  ( YM  -  Yt)  Ym{ }  (2.36) 

where  /  =   1,2  and  XnXM,  Y„  Y,^,  Xm„  Ymt  are  defined  as  before. 
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III.     UNSTEADY  FLOW  PROBLEM  FORMULATION 

A.  GENERAL 

The  modification  work  for  the  unsteady  flow  is  much  more  involved.    It  requires  the 
following: 

1.  The  establishment  of  five  frames  of  reference  viz  one  fixed  inertia  frame  of  reference 
(global),  two  moving  local  frames  of  reference*  and  two  frozen  local  frames  of  ref- 
erenced 

2.  Reformulation  of  the  two  Kutta  conditions  which  are  coupled  non-linearly.  The 
solution  requires  an  iterative  procedure  to  compute  for  the  two  >•(/).  There  are  two 
possible  solutions  to  the  Kutta  condition  due  to  the  quadratic  nature  of  the 
equations.  The  solution  which  ensures  that  the  product  of  the  tangential  velocities 
of  the  first  and  last  panels  of  each  airfoil  is  negative  is  accepted  as  the  solution. 

3.  The  creation  of  a  new  subroutine  (SUBROUTINE  NEWPOS)  which  transforms 
all  coordinates  in  either  of  the  two  respective  local  frames  of  reference  to  the  global 
frame  of  reference.  This  simplifies  the  definition  of  the  airfoil,  wake  element  and 
core  vortices  relative  global  geometries.  The  code  requires  that  the  airfoil  be  de- 
fined with  respect  to  the  respective  frozen  local  frames  of  reference  once  only. 
Subsequent  time  dependent  airfoil  motion,  wake  panel  behaviour  and  core  vortices 
convection  are  computer  generated. 

4.  The  introduction  of  a  more  accurate  method  to  obtain  the  velocity  potential  by 
integrating  the  velocity  over  smaller  panels  on  the  airfoil  without  having  to  store 
large  arrays  of  influence  coefficients  which  are  not  needed  for  satisfying  the  flow 
tangency  condition. 

5.  Extension  of  the  influence  coefficient  to  include  the  effects  of  the  second  airfoil  with 
its  own  peculiar  wake.  This  also  requires  an  introduction  of  an  additional  influence 
coefficient,  that  on  the  wake  element  due  to  the  wake  element  from  the  other 
airfoil. 

B.  FRAMES  OF  REFERENCE 

The  inertia  (global)  and  the  two  frozen  local  frames  of  reference  at  a  specified  time 
tk  are  defined  as  for  the  steady  flow  case.  The  two  moving  local  frames  of  reference  have 
their  x-y  axes  as  for  the  frozen  local  frame  of  reference  but  this  frame  is  moving  with  the 
airfoil. 


4  Moving  local  frames  of  reference  are  used  to  satisfy  the  flow  tangency  equation  which  sim- 
plifies to  equation  2.9. 

5  Frozen  local  frames  of  reference  are  inertia  frames  (used  in  the  steady  flow  case)  of  reference 
used  to  convect  the  core  vortices  of  the  previous  time  steps  with  respect  to  that  time  step  frame  of 
reference  and  subsequently  transformed  to  the  current  time  step  frozen  frame  of  reference. 


C.     UNSTEADY  FLOW  MODEL 
1.     Rigid  Body  Motion 

The  rigid  body  motion  for  the  two  airfoil  system  is  an  extension  to  the  single 
airfoil  system.  Both  airfoils  are  considered  to  have  a  mean  velocity  of—  Vx,  with  a  time 
dependent  translational  velocity  of-[L'(/),/  +  V(f)J]  and  a  rotational  velocity  —  Q(/)* 
which  can  be  in  phase  or  out  of  phase;  i  and  j  are  unit  vectors  in  the  respective  local 
frames  of  reference  and  Q.  is  positive  in  the  clockwise  direction.  The  flow  will  be  deter- 
mined with  respect  to  the  moving  local  frame  of  reference.  The  flow  tangency  condition 
takes  its  simplest  form  in  the  moving  frame  of  reference.  The  flow  tangency  condition 
seen  from  this  frame  of  reference  will  satisfy  equation  2.9.  The  unsteady  stream  velocity, 
Vsrr,.m  is  made  up  by  the  vector  sum  of  a  mean  velocity  Vx,  a  time  dependent 
translational  velocity  [U(l)ki  +  l'(f)J]  and  a  rotational  velocity  Q([)k  (Figure  A). 

V5n«m  =  VJ.  Cos[  y.(f)k  }  7  +  Sin\  *(/)*  ]/}  +  [  U(f)J  +  V{l)J]  +  Q(/)A  (y7  -  xj)        (3.1) 

The  disturbance  potential,  of  necessity,  is  defined  with  respect  to  the  inertia  frame  of 
reference,  for  it  is  only  in  this  frame  that  the  flow  is  irrotational.  It  is  then  transformed 
to  the  moving  frame  of  reference  for  ease  in  treating  the  flow  tangency  condition.  The 
disturbance  potential  is  redefined  to  include  the  contributions  from  the  wake  panel  and 
the  core  vortices  from  both  airfoils. 

6  =  <t>5  +  <t>v  +  4>* .  +  <J)CV  (3.2) 

The  total  velocity  is  then  written  : 

The  unsteady  Bernoulli's  equation  for  the  pressure  coefficients  on  the  airfoil  surface  is 
written  with  respect  to  the  moving  frame  of  reference.  Giesing  showed  this  to  be  written, 
in  our  notation  as  : 

c  ~  (     s'rec"r'  i2  _  i     tolal  i2  _     2     c4)  C3  4) 

"     oc 


oo 


2 
where  Vttnam  and  V,0!a,  are  defined  according  to  equations  3.1  and  3.3  respectively 
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Figure  4.      Frames  of  Reference  for  Airfoil  1 
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2.  Wake  Model 

Recall  that  the  unsteady  flow  model  requires  an  additional  model  to  simulate 
the  continuous  shedding  of  vorticity  into  the  trailing  wake.  The  treatment  of  this  vortex 
shedding  process  follows  the  approach  of  Basu  and  Hancock.  The  shed  vorticity  takes 
place  through  a  small  straight  line  wake  element  attached  as  an  additional  panel  to  the 
trailing  edge  of  each  airfoil  with  a  uniform  vorticity  distribution  [>•„(/)]*•  This  shed 
vorticity  panel  will  be  established  if  its  length  A(/)*  and  inclination  ©(/),  to  the  respective 
local  frames  of  reference  satisfy  the  Helmholtz  theorem  : 

A(/),[y,(0]/.  +  r,(/)  =  r,_1(/)  (3.5) 

or  A  (/)  k  [  yjft  )k  =  r,_,(/)  -  T,(/)  =  SS(f)  [  y(/)A_,  -  y(0*  1  (3-6) 

where  SS  is  the  perimeter  of  the  airfoil  and  I"*.,  and  yk_x  are  respectively  the  total  circu- 
lation and  vorticity  strength  which  are  determined  at  the  previous  time  step  tk_x. 

At  the  next  time  step  r^,,  the  shed  vorticity  panel  will  be  detached  from  the 
trailing  edge  and  will  convect  downstream  as  a  concentrated  free  vortex  with  circulation 
A(/)4  [,'»(/)]i  or  r4_,(/)  —  rk{f)  at  the  resultant  local  velocity  of  the  fluid  particle.  This  wake 
convection  process  is  illustrated  in  Figure  5  where  the  airfoils'  subscripts  are  dropped 
without  loss  of  generality. 

In  the  code,  the  convection  of  the  core  vortices  is  broken  into  three  steps: 

1.  The  core  vortices  are  first  convected  using  the  resultant  absolute  velocity  with  re- 
spect to  the  frozen   local  frame  of  reference  of  the  previous  time  step. 

2.  This  is  followed  by  a  transformation  to  the  current  frozen  local  frame  of  reference. 

3.  Finallv  this  is  transformed  to  the  elobal  frame  of  reference  bv  using  the  new  sub- 
routine (SUBROUTINE  NEVVPOS). 

3.  Additional  Boundary  Conditions 

The  unsteady  flow  model  has  now  introduced  an  additional  boundary  condition 

viz.  the  conservation  of  vorticity  (equation  3.5  or  3.6)  through  the  modeling  of  the  wake. 

However,  the  introduction  of  the  wake  creates  three  additional  unknowns  for  each 

airfoil,  that  is,  yw(f)k,  A(/)*  and  ©(/),.   As  such  two  additional  conditions  are  required  for 

each  airfoil  in  order  to  solve  the  system.  The  approach  suggested  by  Basu  and  Hancock 

is  extended  to  the  two  airfoil  case. 

1.  The  wake  panel  is  oriented  in  the  direction  of  the  local  resultant  velocity  at  the 
panel  midpoint. 
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Figure  5.      Panel  Methods  Representation  for  Unsteady  Flow 


22 


tan  0(/\  =    — : H7) 

where  [VJJf)]«  and  [*£(/)],  are  the  x  and  y  velocity  components  at  the  midpoint  of 
the  wake  panel  with  respect  to  the  frozen  local  frame  of  reference6. 

2.  The  length  of  the  wake  panel  is  propotional  to  the  magnitude  of  the  local  resultant 
velocity  at  the  panel  midpoint  and  the  step  size  of  the  time  step. 

mk  =  &  -  ^-OV'a^oy2 + a^ow2  (3.8) 

We  now  have  all  the  necessary  equations  to  solve  for  the  system. 

D.     INFLUENCE  COEFFICIENTS 

As  with  the  boundary  condition,  additional  influence  coefficients  need  to  be  defined 
as  a  result  of  the  wake  model.  The  definitions  in  the  single  airfoil  case  are  thus  extended 
as  follows: 

1.     More  A's  and  B's  Influence  Coefficients 

Defining  NP3  and  NP4  as  the  panel  number  for  the  wake  for  the  first  and  sec- 
ond airfoil  respectively  and  h  as  a  arbitrary  core  vortex,  the  following  additional  influ- 
ence coefficient  are  required. 

•  (B'yFJ).  :  normal  velocity  component  induced  at  the  i'h  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  wake  panel  of  the  first  airfoil  at  time  tk. 

•  (B ' ::rX  :  tangential  velocity  component  induced  at  the  i'h  panel  control  point  by 
unit  strength  vorticity  distribution  on  the  wake  panel  of  the  first  airfoil  at  time  tk. 

•  {By;,).  :  normal  velocity  component  induced  at  the  /''■  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  wake  panel  of  the  second  airfoil  at  time  tk. 

•  (B:,^).  tangential  velocity  component  induced  at  the  i"1  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  wake  panel  of  the  second  airfoil  at  time  tk. 

•  M'nA  :  x  velocity  component  induced  at  the  first  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  source  distrib- 
ution on  the/"  panel  at  time  tk. 

•  {Afmj)l  :  y  velocity  component  induced  at  the  first  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  source  distrib- 
ution on  the/'  panel  at  time  ih. 

•  (^.'wj*  :  x  velocity  component  induced  at  the  second  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  source  distrib- 
ution on  the/''  panel  at  time  ik. 


6  It  is  important  to  note  that  the  velocities  at  the  midpoint  of  the  wake  panels  do  not  include 
the  effect  due  to  the  vorticity  distributed  along  itself  but  it  does  include  the  contribution  due  to  the 
wake  panel  of  the  other  airfoil. 


23 


• 


• 


• 


• 


(^fcpv)*  '•  >'  velocity  component  induced  at  the  second  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  source  distrib- 
ution on  the/'  panel  at  time  tk. 

{B%-nj)k  '■  x  velocity  component  induced  at  the  first  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  vorticity  dis- 
tribution on  the/'  panel  at  time  tk. 

{Bl-piJ)k  '■  y  velocity  component  induced  at  the  first  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  vorticity  dis- 
tribution on  the/'  panel  at  time  tk. 

(B'SPAJ)k  :  x  velocity  component  induced  at  the  second  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  vorticity  dis- 
tribution on  the/  panel  at  time  tk. 

{&[fp4J)k  :  y  velocity  component  induced  at  the  second  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  vorticity  dis- 
tribution on  the/  panel  at  time  tk. 

{Af,J)k  :  x  velocity  component  induced  at  the  h"'  core  vortex  with  respect  to  the 
frozen  local  frame  of  reference  by  unit  strength  source  distribution  on  the 
/  panel  at  time  ik. 

{A),)\k  :  y  velocity  component  induced  at  the  h"'  core  vortex  with  respect  to  the 
frozen  local  frame  of  reference  by  unit  strength  source  distribution  on  the 
/  panel  at  time  tk. 

(B;;Jk  :  x  velocity  component  induced  at  the  h'h  core  vortex  with  respect  to  the 
frozen  local  frame  of  reference  by  unit  strength  vorticity  distribution  on  the 
/  panel  at  time  tk. 

•  {Bj,)k  :  y  velocity  component  induced  at  the  h'h  core  vortex  with  respect  to  the 
frozen  local  frame  of  reference  by  unit  strength  vorticity  distribution  on  the 
/  panel  at  time  tk. 

•  (Bf,_xn)k  :  x  velocity  component  induced  at  the  h!h  core  vortex  by  unit  strength 
vorticity  distribution  on  the  wake  panel  of  the  first  airfoil  at  time  ik. 

•  {Bi,spi)k  '■  v  velocity  component  induced  at  the  h'h  core  vortex  by  unit  strength 
vorticity  distribution  on  the  wake  panel  of  the  first  airfoil  at  time  ik. 

•  {Bj;WP4)k  :  x  velocity  component  induced  at  the  h'h  core  vortex  by  unit  strength 
vorticity  distribution  on  the  wake  panel  of  the  second  airfoil  at  time  tk. 

•  {Bi^Pi)k  :  y  velocity  component  induced  at  the  h,h  core  vortex  by  unit  strength 
vorticity  distribution  on  the  wake  panel  of  the  second  airfoil  at  time  tk. 

2.     New  C's  Influence  Coefficient 

The  single  airfoil  definition  is  extended  to  include  the  effects  of  the  second 
airfoil. 

•  (Qm(/))*  :  normal  velocity  component  induced  at  the  i'h  panel  control  point  by 
unit  strength  m'h  core  vortex  at  time  th. 
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•  (Q„(/));  :  tangential  velocity  component  induced  at  the  /"'  panel  control  point 
by  unit  strength  m'"  core  vortex  in  the  wake  of  the  /'*  at  time  tk. 

•  (Qra^i(O)*  :  x  velocity  component  with  respect  to  the  frozen  local  frame  of  refer- 
ence induced  at  the  first  airfoil  wake  element  by  unit  strength  m"'  core  vortex  in  the 
wake  of  the  /•'•  at  time  tk. 

•  {QmJJ))k  '■  >' velocity  component  with  respect  to  the  frozen  local  frame  of  refer- 
ence induced  at  the  first  airfoil  wake  element  by  unit  strength  m"'  core  vortex  in  the 
wake  of  the  l'h  at  time  ik. 

•  {Os-Pim(f))k  :  x  velocity  component  with  respect  to  the  frozen  local  frame  of  refer- 
ence induced  at  the  second  airfoil  wake  element  by  unit  strength  m'h  core  vortex  in 
the  wake  of  the  /'"  at  time  tk. 

•  (Qv.F4.m('0)*  :  y  velocity  component  with  respect  to  the  frozen  local  frame  of  refer- 
ence induced  at  the  second  airfoil  wake  element  by  unit  strength  mth  core  vortex  in 
the  wake  of  the  /''  at  time  tk. 

•  (QJJ)),  :  x  velocity  component  with  respect  to  the  frozen  local  frame  of  reference 
induced  at  the  h'n  core  vortex  by  unit  strength  m'h  core  vortex  in  the  wake  of  the  /"' 
at  time  /,. 

•  (Q  J/));  :  y  velocity  component  with  respect  to  the  frozen  local  frame  of  reference 
induced  at  the  h*  core  vortex  by  unit  strength  m"'  core  vortex  in  the  wake  of  the  l"1 
at  time  i>. 

3.     Computation  of  the  Time-Dependent  Influence  Coefficient 

As  for  the  single  airfoil  ,  the  influence  coefficients  of  the  wake  element 
(B\\n)*'{Blvn)K(B'.yp-i)  an^  (^v/*)*-  are  computed  as  for  the  airfoil  influence  coefficients 
B"  and  B'  with  the  subscripts  NP3.NP-4  replacing  j. 

The  x  and  y  velocity  components  for  the  respective  frozen  local  frames  of  ref- 
erence are  obtained  indirectly  by  first  calculating  for  the  global  X  and  Y  velocity.  These 
global  velocities  are  then  transformed  to  the  frozen  frames  through  a  simple  relationship 
by  the  use  of  the  respective  angle  of  attack.  It  can  be  easily  shown  that  the  velocity  with 
respect  to  the  frozen  frame  of  reference  is  related  to  the  global  velocity  by  the  following 
relationship. 

Kmh  =  [GV-m  cos  <x(l)]k-[GVlm  sin  «(/)], 

(Km)k  =  [GV*msma([)}k  +  {GVlmcosa.([))k  (3-9) 

where  V  denotes  a  generic  induced  velocity  with  respect  to  the  frozen  frame  and  the 
precedent  G  denotes  global  velocities. 
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The  global  {GA'NnJ)k,{GAi-n])k,{GA'NP^)k, ...  {GRyAJ)k,{GA$k,{GA%{GB$k  and 
(Gfyj),,  are  computed  using  equations  2.18  through  2.21  with  0,  set  to  zero  and  subscript 
i  appropriately  replaced. 

The  global  C's  coefficients  will  be  computed  identically  as  for  the  single  airfoil 
case.   The  results  are  repeated  here  for  the  sake  of  completeness. 

{GCim{l))k  =    -  (3.10) 

Z7tVim)k 

.                        Sm\mk  -  (0m)k] 
(GCrim(f))k=    ~   \"  (3-11) 

where  : 


(O,  =  s/Kxm,  -  xj-  +  [y-m,  -  ymf\ 

xm,  =  1/2  (x,  +  x„x) 

ym,=  1/2  Cv,  +>•_,) 

xm  =    x  coordinate  of  m"'  core  vortex  at  time  ik 

y„  =    y  coordinate  of  m"'  core  vortex  at  time  tk 

0,  =  ian~\  }x'lZ}x   ) 

■*i-]        -s 

vm,  —  v_ 

Also    (GC^3,m(/))*,  (GCW(0)*.(GC3^(0)to(GQ^(/))*,(GC5U())*    and    (GQ,m(/))<    are 

computed  with  equation  3.10  with  6,  set  to  zero  and  the  subscript  i  appropriately  re- 
placed. 

E.     NUMERICAL  SOLUTION  SCHEME 
1.     The  Flow  Tangency  Condition 

The  unsteady  flow  tangency  equation  for  the  two  airfoil  system  is  a  simple  ex- 
tension of  the  single  airfoil  case. 

n  In  n  In 

Y}AfaM + Yj  [A'j{qj)k] + y{i)kYjB*j + y{2)k  Yj  b» + [( ^siream)i '  "i]k 

j=\  j=n+\  j=l  j=n+] 

+  [>'J1)]/Xyp3),<  +  [Y&MBZnJt  + 
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k-\  k-\ 


U:=l  J  Ui=l  J 

where  i=  1.2  ...  2n 

(3.12) 

where  i=  1.2  ...  2n  and  Vstrtam  is  evaluated  by  equation  3.1  at  the  /"'  panel  control  point. 

Rearranging  the  equation  we  can  express  the  source  strengths  explicitly  as  a 
function  of  the  two  vorticity  strengths  and  a  constant  term  as  in  equation  2.28  where  : 


n  2n 

L.H.S.  =  2jAfaj)k]+  2ji*jfa)jkJ  <3'13) 

n 

First  R.H.S.  =  y(l)J(  -^yp  )k{B^P3)k  -  ]m/?J)A  (3.14) 

S^nrf  R./f.S.  =  y(2)  J(  -^-  )a(2£vM)a  -  ^  (J^fcj  (3.15) 

-*  SS(1) 

ThndR.H.S.   =     -[(Vstream)l'ni\k-yV)k-l(-^MB?,SP3)k- 


k-\  k-\ 


2j{c"Mwm-M  -  rm(/))]L,  -  2,^c^)(r-»(/) "  r«W)l[w  (3-i6) 

where  i  =   1,2  ...  2n 

For  the  code,  the  above  is  implemented  in  a  2A'x(2Ar+3)  matrix  in  SUB- 
ROUTINE COEF  and  is  subsequently  solved  by  GAUSS  elimination  to  obtain  an  ex- 
plicit expression  for  the  source  strength  in  terms  of  the  vorticity  distribution  and  a 
constant  term  as  in  the  steadv  flow  case  as  follows  : 
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qj  =  b\j.y(l)  +  b2j.y{2)  +  b3j   ,  j  =  1,2 ...In  (3.17) 

where 


part  of  the  source  strength  which  depends  on  circulation  y(l)A. 
part  of  the  source  strength  which  depends  on  circulation  y{2)k. 
part  of  the  source  strength  which  is  independent  of  the  circulation. 


Equation  3.17  is  critical  in  the  simplification  of  the  treatment  for  the 
non-linear  coupled  Kutta  condition. 

2.      The  Kutta  Condition 
For  the  unsteady  flow,  the  Kutta  condition  can  be  written  as  the  following 
For  the  first  airfoil 

(K)*-('»)*  =  2[  47(4>n-4>,))k  =  21-^T1)* 


h  ~  h-i 


For  the  second  airfoil 

-255(2)  ^l;'-^ 


h  -  lk-\ 


The  tangential  velocity  for  the  first  panel  can  be  written  thus: 


In 


(3-18) 


(3-19) 


v\  =  -  Y}Bl^  +  yWkYjAv +  +  y{2)k  Yj  A"j  + 

;=i  j=\  J=n+\ 

•^y-  <Ar/.3[y*-iO)  -  y±(i)]  +  -^-  ^UJy*-i(2)  -  y*(2)]  +  K^J,  •  ?]*  + 

ic-l  A-l 

\Yfl<ZMXrm-l®  ~  rmW)]}w  +  JSK^-W^-lW  ~  r«WMJw  (3-20) 

l»i«i  J         Ui=i  J 

Substituting  equation  3.17  into  3.20,  the  following  simplified  relation  for  V\  is  obtained. 
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V[  =  D\,y(l)k  +  D2xy{2)k  +  Dl,  (3.21) 

where 


2n  n 


02,=  -][^2l+£4_^|L<w 


7=1  y'=n+l 


ESSiX)  SS(~>)  -  - 

7=1 
*-l  Jt-1 

+  jyjfcUoxr,^/)  -  rm(/))]j/«,  +  j^Kc^oxr-,.^  -  rj/))l}/=2  (3-22) 

Similar  expressions  can  be  obtained  for  the  other  trailing  edge  panels  as  follows  : 

r-  =  Dlny{\)k  +  D2ny{2)k  +  D3n  (3.23) 

K+i  =  ^In+iKl)*  +  D2n+ly(2)k  +  D3M+1  (3.24) 

V'2n  =  D\2ny(\)k+D22nV(2)k  +  Dhn  (3-25) 

where  the  coefTicients  Dl,  D2,  D3  are  obtained  as  per  equation  3.22.    We  require  the 
square  of  the  tangential  velocity  for  the  Kutta  condition 

(I'tf  =  (D\,y(l)k  +  D2,y(2)k  +  D3,)2 

=  Dl?y(l)J  +  D2?y(2)J+Z)3;  (3.26) 

+  2Z>l1Z)2Iy(l)*y(2)i  +  2/>l1Z)31y(l)i  +  2Z>21Z)31y(2), 

(V'n)2  =  (D\ny(\)k  +  D2ny{2)k  +  D3n)2 

=  Dl2ny{\)2k  +  D22ny{2Jk  +  D32n  (3.27) 

+  2DlnD2„y(l)Ay(2)A  +  2D\nD3ny{\)k  +  2D2nD\y{2)k 
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For  the  first  airfoil  equations  3.26  and  3.27  are  then  substituted  into  3.18  to 
give,  for  the  code: 

AAA(\)y(\)2k  +  BBB(\)y{2)2k  +  CCC(l)y(l)k  + 

DDD(l)y(2)k  +  EEE(\)y(l)ky(2)k  +  FFF(l)      =    0  (3.28) 

where 

AAA{\)  =  D\]-D\2n 

BBB{\)  =  Dl]-D22n 

SS(\) 
CCC(l)  =  2[Dll.D3,-D\n.D3n-      _)>    ] 

lk       lk-\ 

DDD(l)  =  2[D2,  .  Dlx  -  D2n .  Z)3„] 
EEE{\)  =  2[D\l*D2x-Dln»D2n] 

2  2  SS(1>A_,(1) 

FFF(\)  =  Dl\  -D?n  +  2(  k  )  (3.29) 

'k       [k-\ 

A  similar  set  of  equation  can  be  obtained  for  the  second  airfoil  as  follows: 
AAA(2)y(\)2k  +  BBB{2)y(2)l  +  CCC(2)y(\)k  + 

DDD(2)y(2)k  +  EEE(2)y(\)ky(2)k  +  FFF(2)      =    0  (3.30) 

where 
AAA{2)  =  Dl2n+1-D\\n 

BBB(2)  =  D22n+,  -  D2\n 

CCC(2)  =  2[Dln+1.D3n+1-Dl2n.D32n] 

DDD(2)  =  2[D2n+]  .  D3n+]  -  D2ln  .  D32n  -     5f(2)     ] 

lk       'k-\ 

EEE(2)  =  2[Dln+l.D2n+l-D\2n.D22n) 
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FFF(2)  =  Z)3;+1-Z)3^  +  2(5^(2)/;:-l(2)  )  (3.31) 

The  above  Kutta  equations  are  non-linear  and  coupled.  To  linearize  it,  we  adopt 
Newton's  method,  in  which  the  initial  iteration  is  linearized  about  the  values  of  the 
previous  time  step: 

KDS  =  KDCJ  +  Wfi*WiSf  =  [y(i)*-if  +  2y(ifi:i*y(i)2 

y(2)J  =  y{2)T_\  +  5y(2)J  -  [y(2gf  =  [y(2)J_1]2  +  2y(2)£|<5y(2)J  (3.32) 

where: 

^/(D^y(i)r1  ;  Ki)J  =  y(i)*-, 

5y(2)J«y(2)J-1    I    y{2)l  =  y{2)k_, 

with  pc  denoting  the  iteration  counter.  After  dropping  quadratic  terms  of  Sy  and  col- 
lecting like  terms  the  Kutta  condition  for  the  first  airfoil  simplifies  to  the  following: 

{2AAA(l)y(l)k_^  CCC(l)  +  EEE(l)y(2)k_]}Sy(l)k  + 
{2BBB(\)l(2)k_,  +  DDD{\)^-  EEE(\)y(\)k_l}Sy(2)k  + 
AAA[\)y{\)\_,  +  BBB{\y,i2)l_,  +  CCC(1)>-(1)A_,  + 
DDD{\)y{2)k_^  EEE[\y;{\)K_,y{2)k_,+  FFF{\)  =  0 

(3.33) 
Similarly  the  linearised  Kutta  condition  for  the  second  airfoil  can  be  written: 
{2AAA(2)y(l)k_l  +  CCC(2)  +  EEE(2)y(2)k_]}Sy(l)k  + 

{2BBB(2)y(2)k_,  +  DDD(2)  +  EEE(2)y(l)k_l}Sy(2)k  + 
AAA{2)y{\)\_,  +  BBB(2)y(2)2k_]  +  CCC(2)y(l)^,  + 
DDD(2)y(2)k_l  +  EEE(2)y(\)k_ly(2)k_,+FFF(2)  =  0 

(3.34) 
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The  above  two  linear  equations  with  two  unknowns  Sy(l)k  and  6y{2)k  can  now  be  easily 
solved  with  the  same  Gauss  routine.  The  results  are  then  back  substituted  into 
equations  3.32.  This  procedure  is  repeated  until  the  corrections  dy(l)j  and  6y(2)*k  are  less 
than  a  prescribed  tolerance. 

3.  Convection  of  the  Wake 

The  resultant  velocities  of  all  core  vortices  are  calculated  using  the  frozen  local 
frame  of  reference.   This  is  resolved  into  components  at  the  h"'  core  vortex  as  follows: 

n  2n  n  2n 

(Uh)k  =  2JAhMA  +  Yj  [^)]*  +  vVh^JBtik  +  y(2)k  ^  (Bxhj)k  +  [{VJh  .7}k 

j=\  j=n+]  J=\  J=n+\ 

fc-i  fc-i 

+  ly.KcUour^if)  -  rm(/))]j/=1  +  J^KC^^r^/)  -  rjOfljw 

m^h  m=th 

(3.35) 

n  In  n  In 

{Vh)k  =  2JAyhj{qj)]k  +  ^T  [Ayhj{qj))k  +  y(l);^(4),  +  y(2),]T {Byj)k  +  [(VJh  .j]k 

J=]  ;=«+]  j=\  j-n+l 

+  [vM]k(BLNp3)k  +  [y*(2)]k(BlSPA)k 

k-\  k-\ 

+  ^[(cUOU^nUO  -  rm(/))lL,  +  jX[(cL(0)*(rm_,(/)  -  rm(/))]L2 

m^h  m&h 

(3.36) 

The  location  of  the  core  vortices  at  the  new  time  step  is  then  computed  with  respect  to 
the  current  time  step  frozen  local  frame  of  reference  and  then  transformed  to  the  new 
time  step  frozen  local  frame  of  reference.  These  coordinates  are  subsequently  trans- 
formed to  global  coordinates  so  as  to  facilitate  the  calculation  of  the  influence  coeffi- 
cient. 

4.  Disturbance  Potential  and  Pressure  Distribution 

The  concept  of  disturbance  potential  introduced  in  the  steady  flow  has  played 
a  vital  role  in  that  it  facilitates  the  synthesis  of  the  flow  field  from  simple  flow  field  by- 
simple  superposition  of  the  various  singularities  contributions.  However,  there  has  not 
been  a  requirement  to  solve  directly  for  the  disturbance  potential,  as  our  interest  was 
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on  the  disturbance  induced  velocity  which  is  the  spatial  derivative  of  the  disturbance 
potential.  In  our  solutions,  the  concept  of  influence  coefficient  has  allowed  a  direct 
evaluation  of  the  disturbance  velocity  thus  nullifying  the  requirements  of  obtaining  the 
disturbance  potential. 

The  treatment  for  the  unsteady  flow,  though  it  follows  the  same  procedure  as 
the  steady  flow  case,  viz  the  influence  coefficient  to  obtain  the  disturbance  velocity,  still 
requires  the  disturbance  potential  for  the  computation  of  the  pressure  coefficient  as  can 
be  seen  in  equation  3.4.    Rewriting,  we  have: 

\<-p)k  =  l(  ^        h\k  ~  K  HP —  hik  -~2 7-r~, 1-1,2  ...  In  (3.37) 


where  \\„tcm  and  Vlotal  are  defined  according  to  equations  3.1  and  3.3  respectively. 

From  equation  3.37  we  have  written  the  rate  of  change  of  4>  by  a  backward  fi- 
nite difference  approximation.  This  simplifies  our  iteration  procedures  tremendously  as 
the  <t>  from  the  previous  time  step  exists  at  the  current  computation.  The  computation 
of  the  disturbance  potenial  <b  is  obtained  through  two  steps.  The  difference  in  potential 
<t>  from  upstream  at  infinity  to  the  leading  edge  is  computed  and  is  then  combined  with 
the  difference  in  the  potential  from  the  leading  edge  to  the  panel's  of  interest  control 
point. 

The  present  approach  differs  from  the  original  single  airfoil  approach  in  that  the 
velocity  potential  along  the  airfoil  surface  is  computed  via  a  finer  grid  using  Gaussian 
quadrature".  This  modification  improves  the  results8  but  the  cost  of  doing  it  is  a  longer 
computation  time.  The  definition  of  infinity  has  also  been  set  to  100  chord  lengths  in 
comparison  to  10  chord  lengths  used  in  the  original  code.  However,  the  number  of 
computation  points  and  the  first  panel  length  from  the  leading  edge  upstream  have  not 
been  changed. 

For  completeness,  the  computation  of  the  disturbance  potential  from  the  lead- 
ing edge  to  infinity  is  included  here  for  both  airfoils.  It  is  essentially  the  same  as  the 
single  airfoil  case  except  that  the  disturbance  potential  is  now  computed  relative  to  the 
global  frame  of  reference  instead  of  the  airfoil  fixed  frame  of  reference  as  used  in  the 


7  Each  panel  is  subdivided  into  4  additional  sub-panels  and  the  tangential  velocity  is  computed 
and  integated  over  these  smaller  panels  with  a  weighting  function  to  get  the  disturbance  potential. 

8  One  can  check  the  improvement  by  computing  the  y(s)  obtained  through  the  difference  be- 
tween the  trailing  edges  and  compares  them  with  that  obtained  through  the  Kutta  condition. 
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original  code.    We  begin  by  selecting  a  straight  line  extending  upstream  in  the  direction 

parallel  to  Vx.  The  length  of  the  line  is  set  at  100  chord  lengths.  This  line  is  divided  into 

z  panels  with  the  first  element  at  the  leading  edge  set  equal  to  the  single  airfoil  case  for 

the  purpose  of  ensuring  that  the  panel  size  is  comparable  to  the  airfoil  panel  size.   The 

panel  size  is  then  subsequently  increased  to  take  advantage  of  the  inversely  decaying 

induced  velocities  at  larger  distances.  Using  subscript  f  to  denote  these  panel  mid-points, 

we  define  the  following  influence  coefficients: 

{ArfJ)k    :  normal  velocity  component  induced  at  the  fh  panel  control  point  by  unit 
strength  source  distribution  on  the /''panel  at  time  tk. 

{A'A,)k  :  tangential  velocity  component  induced  at  the/"1  panel  control  point  by  unit 
strength  source  distribution  on  the  /''panel  at  time  tk. 

(#;,)*    :  normal  velocity  component  induced  at  the/''  panel  control  point  by  unit 
strength  vorticity  distribution  on  the /''panel  at  time  tk. 

{B}x)k  :  tangential  velocity  component  at  the/''  panel  control  point  by  unit  strength 
vorticity  distribution  on  the /''panel  at  time  tk. 

(Bf,.\'pi)k  '■  normal  velocity  component  induced  at  the/''  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  wake  panel  of  the  first  airfoil  at  time  tk. 

(Bfjrn)*    '■    tangential  velocity  component  induced  at  the/'1  panel  control  point  by 
unit  strength  vorticity  distribution  on  the  wake  panel  of  the  first  airfoil  at  time  tk. 

(B/.spJk  '•  normal  velocity  component  induced  at  the/''  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  wake  panel  of  the  second  airfoil  at  time  tk. 

{B'fsSP^k  tangential  velocity  component  induced  at  the/''  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  wake  panel  of  the  second  airfoil  at  time  tk. 

(Qm(/))A    :  normal  velocity  component  induced  at  the/''  panel  control  point  of  the 
/'''  by  unit  strength  mth  core  vortex  at  time  tk. 

{Cym{t))k     :    tangential  velocity  component  induced  at  the/'1  panel  control  point 
of  the  /■''  by  unit  strength  m'h  core  vortex  at  time  ik. 

The  tangential  velocity  at  the/''  panel  written  as  for  the  code  is  then: 


in  n  in 

(*/)*  =  -  Yj{B^j)k + y(i)*X(/^* + + y(2)*  Yj  {A*)k + 

j=\  ;=i  j=n+\ 

~^JJ~  (^vP3)fc[^-i(l)  -  V*(l)]  +    A/2)    (4W*[v*-i(2)  -  yk{2)}  + 

k-\  k-\ 

{Xl(c/m(0)A(rm_1(0  -  rm(/))]L  +  (X[(c;m(/))*(r,*_1(')  -  rm(0)]jw  (3-38) 


<m=\  '  K-m=\ 
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valid  for  f=  1,2  ...  z.    The  disturbance  potential  at  the  airfoil  leading  edge  is  the  sum  of 
the  products  of  the  disturbance  induced  velocity  at  each  panel  and  the  panel  length  . 

z 

(W*  =  -  Z[( l%  •  ^  -  xfi  (3-39) 

/=1 

The  integral  over  the  airfoil  surface,  as  stated  before,  is  now  done  over  a  finer 
grid.  This  requires  the  computation  of  the  disturbance  induced  velocity  over  the  smaller 
grids  within  each  panels.  Defining  the  total  finer  grid  points  in  one  panel  by  P9,  we  de- 
fine first  the  refined  influence  coefficient 

tangential  velocity  induced  at  the  /'*  panel  p'h  node  due  to  unit  strength 
source  distribution  on  the/*  panel  at  time  tk 

The  other  refined  influence  coefficients  have  the  same  definition.    The  tangential  com- 
ponent of  the  disturbance  induced  velocity  at  the  i'h  panel  m'h  node  is  then 


In 


(V'ip)k  =    -^(^),  +  y(l)^^  +  y(2)^^  + 

7=1  7=1  j=n+\ 


*-i  *-i 


(OXr^O  -  rw(/))]L,  +  {XkqWWX^CO  -  rm(0)]L2  (3.40) 


^r>i=]  ^  ^m=) 


which  is  valid  for  i=  1.2  ...  2n  ;  p=  1,2  ...  P.    Performing  the  line  integration  by  summa- 
tion, the  disturbance  potential  at  the  i*  nodal  point  is  then 

Wnode  l)k  -  (<*>le)*  +    )    Xl(  VlJ*JJ+l  Wp\  »     /*'  w  >  '  ^  'te 

4-'/?=i 

(3.41) 
'^j  p 

=  (0|e)*  ~    /    Z[( ^Wfl/H  "P]  '     f°rn>i>  <le 


W/*1 


9  This  includes  the  end  points. 
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where  rJJJ.x  denotes  the  panel  length  and  W  the  Gaussian  quadrature  weighting  function. 


rJJ+]  =  v (A}+1  -  X/  +  (}}+1  -  }})2  (3.42) 

Finally  the  disturbance  potential  at  the  i'h  panel  control  point  is  : 

M)*   =    U2[{4>nodel)k  +  (<t>nodei+l)k]    >    <   =    U ...  2«  (3.43) 

5.     Computation  of  Forces  and  Moments 

The  C/,  Cd  and  Cm  about  the  leading  edge  of  each  airfoil  are  calculated  in  exactly 
the  same  way  as  it  is  done  for  the  steady  flow  problem  by  integrating  the  pressure  dis- 
tribution (See  section  D-4  of  Chapter  2). 
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IV.     DESCRIPTION  OF  COMPUTER  CODE  USPOTF2 

A.     PROGRAM  STRUCTURES  AND  CAPABILITIES 

1.  Restrictions  and  Limitations 

The  restrictions  and  limitations  listed  in  the  single  airfoil  documentation  [Ref. 
1]  still  apply  for  the  two  airfoil  case.  Some  efforts  were  made  to  optimise  storage  re- 
quirements versus  computational  repetitions;  however  as  in  the  original  code  there  is  still 
much  room  for  improvement.  Krainer  [Ref.  6  ]  has  improved  the  original  code  by 
combining  subroutines  and  reducing  repeated  computation  by  introducing  additional 
common  variables.    Some  of  his  improvements  have  been  added  to  this  program. 

The  computer  system  used  is  the  Naval  Postgraduate  School  IBM  3033AP.  The 
current  program  fixes  the  maximum  number  of  airfoil  panels  to  20010  and  the  maximum 
allowable  time  steps  to  200.  The  computer  currently  has  to  be  run  with  a  minimum 
storage  requirement  of  2  Mbytes.  A  detailed  computing  time  study  for  the  program  is 
not  undertaken,  for  it  not  only  changes  with  the  number  of  nodal  points  selected  but 
also  with  the  time  step  increment  as  this  has  an  effect  on  the  rate  of  convergence  for  the 
wake  panel  iteration.  An  order  of  magnitude  is  given  for  one  case  run  to  give  an  ap- 
preciation of  the  computing  time  required.  The  system  currently  requires  a  total  CPU 
run  time  of  200  seconds  using  an  optimising  compiler  for  a  step  input  with  a  0.025 
non-dimensionalised  time  increment  for  26  time  steps. 

At  present,  the  program  can  run  for  two  airfoils  set  at  arbitrary  distance  and  at 
different  angles  of  attack  undergoing  any  of  the  following  motions: 

1.  In-phase  and  out-of-phase  Step  Input 

2.  In-phase  and  out-of-phase  Modified  Ramp  Input 

3.  In-phase  and  out-of-phase  Translational  Harmonic  Oscillation 

4.  In-phase  and  out-of-phase  Rotational  Harmonic  Oscillation 

5.  Sharp  Edge  Gust  Field  Penetration 

2.  Current  Structures  of  USPOTF2  Main  Program 

The  flow  logic  for  USPOTF2  (Unsteady  Potential  Flow  for  2  airfoils)  is  illus- 
trated in  Figure  6.    It  is  essentially  divided  into  three  major  modules,  namely,  the  input 


10  These  are  combined  total  panels  for  both  airfoils. 
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-  problem  setup,  Steady  flow  solution  with  its  associated  output  and  the  Unsteady  flow 
solution  with  its  associated  output. 

The  first  module  includes  subroutines  IXDATA,  SETUP,  BODY  and  NACA45. 
The  main  functions  of  this  module  are: 

1.  Set  up  the  problem  formulation  by  reading  in  the  necessary  values  and  flag  setting 
from  filecode  1. 

2.  If  the  flag  for  a  NACA  4-digit  or  5-digit  of  type  230XX  is  set,  this  module  calls 
subroutines  BODY  and  NACA45  to  compute  the  local  panel  coordinates. 

3.  If  the  flag  for  a  NACA  4-digit  or  5-digit  of  type  230XX  is  not  set,  this  module  then 
proceeds  to  read  in  the  local  panel  coordinates  from  filecode  1. 

4.  The  local  slopes  and  the  airfoil  perimeters  are  then  calculated  in  preparation  for  the 
next  module. 

The  second  module  calculates  the  Steady  flow  solution.  It  calls  subroutines 
NEWPOS.  INFL,  COEF,  GAUSS,  KUTTA,  VELDIS  and  FANDM.  The  main  func- 
tions of  this  module  are: 

1.  Transform  the  local  airfoil  coordinates  into  global  coordinates  and  compute  the 
influence  coefficients. 

2.  Set  up  flow  tangency  equation  as  per  equation  2.27  and  solve  for  the  source 
strengths  (q.)  as  a  function  of  the  vorticitv  strengths  [>*(/)]  and  a  constant  part  where 
j  =  f,2  ...  2nand/=  1,2. 

3.  Set  up  and  solve  the  Kutta  condition  as  per  equations  2.25  and  2.26  for  the 
vorticitv  strengths  and  back  substitute  to  get  the  source  strengths. 

4.  Compute  the  total  tangential  velocity  in  the  moving  frame  and  compute  the  dis- 
turbance potential  (</>,)  and  pressure  coefficient  [(Cp),]  (i  =  1.2  ...  2n). 

5.  Finally,    compute    the    aerodynamic    coefficients    of   forces    and    moments    - 

c  c   c 

6.  This  program  can  terminate  in  this  module  after  all  the  steady  flow  parameters  are 
obtained  without  necessarilv  running  the  Unsteadv  flow  code. 

The  third  module  calculates  the  Unsteady  flow  solution.  It  calls  subroutines 
NEWPOS,  INFL,  COEF,  GAUSS,  KUTTA,  TEWAK,  PRESS,  FANDM  and 
CORVOR.   The  main  functions  of  this  module  are: 

1.  For  the  particular  unsteady  motion,  to  compute  the  initial  time  step  and  new  airfoil 
orientation  with  its  associated  local  body  velocities. 

2.  Introduce  the  wake  panel  and  assume  an  initial  length  and  orientation  to  begin  it- 
erations. 

3.  Transform  all  local  coordinates  and  panel  slopes  to  global  coordinates  and  global 
panel  slopes. 
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4.  Update  influence  coefficient  and  set  up  flow  tangency  equation. 

5.  Solve  for  the  source  strengths  in  terms  of  the  vorticity  strengths  and  the  constant 
part. 

6.  Invoke  the  required  Kutta  condition  and  solve  for  the  vorticity  strengthsil;  back 
substitute  to  get  the  source  strengths. 

7.  Update  the  local  velocities  of  the  wake  element  and  compute  new  length  and  ori- 
entation —  check  for  convergence. 

8.  If  wake  element  velocities  have  not  converged,  iterate  with  the  new  wake  element 
geometries  until  convergence! 2. 

9.  Compute  the  total  velocity,  disturbance  potential  and  the  pressure  coefficient. 

10.  Compute  the  aerodynamic  lift,  drag  and  moment  coefficient. 

11.  Adjust  time  step  either  through  an  interative  inputs  or  by  a  automatic  time  in- 
crement. 

12.  Compute  resultant  local  velocities  of  the  core  vortices. 

13.  Convect  the  core  vortices  by  the  procedure  described  in  Chapter  III  section  E-3 


1 1  There  are  two  possible  solutions;  only  the  vorticity  distribution  that  ensures  that  the  prod- 
uct of  the  relative  tangential  velocities  of  the  upper  and  lower  panels  of  the  trailing  edge  is  negative 
is  accepted  as  solution. 

12  The  convergence  criterion  is  user  specified  through  input  data  for  TOL. 

13  This  is  accomplished  by  setting  TADJ  to  be  non-zero  and  is  intended  for  use  in  conjunction 
with  the  viscous  flow  program. 


39 


READ  INPUT  FROM 
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Figure  6.      Flow  Chart  for  USPOTF2  Computer  Code. 


40 


TRANSFORM  TO  GLOBAL 
COORDINATES  &  SLOPES  & 
COMPUTE  INFL  COEF 

SET  FLOW  TANGENCY  EQN  AS 
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AS  PER  EQN  2.  27 
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SOLVE  FOR  SOURCE  STRENGTH 
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STEADY  FLOW  MODULE 


Figure  6  (Cont'd) 
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INITIALISE  PARAMETERS  BEFORE 
STARTING  UNSTEADY  FLOW  SOLUTION 
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Figure  6  (Cont'd) 


42 


B.     DESCRIPTION  OF  SUBROUTINES 

1.  Subroutine  BODY 

This  subroutine  is  called  for  the  purpose  of  obtaining  the  local  (x,y)  coordinates 
of  either  a  NACA  XXXX  or  230XX  type  airfoil.  It  is  called  by  subroutine  SETUP  and 
it  in  turn  calls  subroutine  XACA45  to  obtain  the  airfoil  thickness  and  camber  distrib- 
utions. 

2.  Subroutine  COEF 

This  subroutine  was  originally  intended  for  the  unsteady  flow  calculation  for  the 
single  airfoil  case.  It  is  now  modified  to  include  the  steady  flow  calculations.  Its  pur- 
pose is  to  set  up  the  flow  tangency  matrix  as  in  equation  2.27  for  the  steady  flow  and 
equations  3.13  through  3.16  for  the  unsteady  flow  case.  These  matrices  are  necessarily 
set  up  in  this  way  so  that  the  source  strengths  can  be  solved  in  terms  of  the  vorticity 
strengths  and  a  constant  by  subroutine  GAUSS  as  a  linear  system  with  three  right  hand 
sides.    It  is  called  by  the  MAIN  program. 

3.  Subroutine  COFISH  (deleted  for  the  t>vo  airfoil  case) 

This  subroutine  was  originally  set  up  to  serve  the  same  function  as  subroutine 
COEF  for  the  steady  case.    It  is  now  deleted  for  computational  efficiency. 

4.  Subroutine  CORVOR 

This  subroutine  is  called  for  the  purpose  of  obtaining  the  convective  velocities 
for  all  the  wake  core  vortices  with  respect  to  the  frozen  local  frame  of  reference  of  the 
current  time  step  in  accordance  with  equations  3.35  and  3.36  where  all  the  influence 
coefficients  are  now  locally  computed  (previously  done  in  subroutine  INFL)  to  save 
storage  requirements  since  these  are  only  required  in  this  subroutine.  The  local  influence 
coefficients  are  obtained  indirectly  by  first  obtaining  the  global  influence  coefficients  and 
transforming  them  viz  equation  3.9.  This  subroutine  is  called  by  the  main  program 
nearing  the  end  of  the  unsteady  flow  calculations  before  starting  a  new  time  step. 

5.  Subroutine  FANDM 

This  subroutine  is  intended  to  calculate  the  overrall  lift  coefficient,  drag  coeffi- 
cient and  moment  coefficient  about  the  leading  edge  for  the  two  airfoils.  In  order  to 
preserve  the  option  of  obtaining  the  x  and  y  forces  with  respect  to  the  respective  local 
frames  of  reference,  the  original  method  rather  than  equations  2.34  through  2.36  is  im- 
plemented. This  subroutine  is  called  by  the  MAIN  program  in  both  the  steady  and  un- 
steady flow  computations  immediately  after  computing  the  pressure  coefficient. 


43 


6.  Subroutine  GAUSS 

This  subroutine  is  the  standard  linear  system  solver  that  employs  the  well- 
known  Gaussian  elimination  with  partial  pivoting  and  operates  simultaneously  on  a  user 
specified  number  of  right-hand-sides.  It  is  called  by  the  MAIN  program  in  both  the 
steady  and  unsteady  flow  calculations.  In  order  to  use  GAUSS,  the  coefficients  of  the 
augmented  matrix  must  be  set  up  so  that  GAUSS  will  return  the  solutions  replacing  the 
corresponding  columns  of  the  augmented  matrix  that  were  initially  occupied  by  the 
right-hand-sides.  The  coefficient  set-ups  are  done  by  subroutine  COEF  for  both  steady 
and  unsteady  flow  problems. 

7.  Subroutine  INDATA 

This  subroutine  is  intended  to  read  in  the  first  three  to  five  cards  of  the  input 
data  depending  on  whether  I  FLAG  =  0.  The  first  three  cards  contain  some  description 
of  the  airfoil  type,  problem  definition,  I  FLAG  information  as  well  as  the  number  of 
lower  and  upper  panels.  If  I  FLAG  =  0,  it  will  treat  the  airfoils  as  NACA  type  airfoils 
and  will  proceed  to  read  the  NACA  number  and  to  calculate  the  thickness  parameters 
that  will  be  required  by  subroutine  NACA45.  This  is  the  first  subroutine  called  by  the 
MAIN  program. 

8.  Subroutine  INFL 

This  subroutine  generates  most  of  the  influence  coefficients  that  are  needed  and 
shared  by  the  different  subroutines.  It  has  been  modified  to  include  the  steady  flow  case 
for  the  purpose  of  reducing  repeated  computations.  It  utilises  the  known  relative  ge- 
ometrical parameters  of  the  singularities  to  carry  out  computation  based  on  equations 
2.18  through  2.21  for  the  steady  flow  calculations  and  including  3.9  through  3.11  for  the 
unsteady  flow  case.  The  MAIN  program  calls  this  subroutine  in  every  iteration  cycle 
of  each  time  step  so  that  the  time  dependent  influence  coefficients  can  be  updated  as 
and  when  neccessary.  Time  independent  coefficients  are  computed  once  in  the  entire 
flow  solutions!-*. There  are  also  time  dependent  influence  coefficients  that  are  independ- 
ent of  the  iterative  cycle  to  obtain  the  wake  panel  orientation^.  Those  influence  coef- 
ficients involving  the  wake  core  vortices  are  also  independent  of  the  iteration  cycle  and 
are  also  updated  once  in  each  time  step;  but  the  process  of  the  update  is  more  compli- 
cated. This  is  due  to  the  process  that  the  wake  need  first  to  be  convected  with  respect 
to  the  previous  time  step  frozen  frame  of  reference,  transformed  to  the  present  frame  of 


14  These  are:  on  the  airfoil  panel  by  the  panels  on  the  same  airfoil. 

15  These  are:   on  the  airfoil  panel  by  the  panels  on  the  other  airfoil. 
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reference  and  finally  transformed  to  the  global  frame  of  reference  by  subroutine 
NEWPOS.  Finally  the  influence  coefficients  involving  the  wake  panels  are  calculated 
as  frequently  as  the  number  of  iterations  take  to  find  a  converged  solution. 

9.  Subroutine  KUTTA 

This  subroutine  is  intended  to  solve  the  Kutta  equation.  It  has  been  modified 
to  include  the  steady  flow  solution.  After  the  source  strengths  have  been  determined  by- 
subroutines  COEF  and  GAUSS  in  terms  of  the  vorticity  strengths  and  a  constant,  this 
subroutine  invokes  the  Kutta  condition  as  in  equations  2.25,  2.26  for  the  steady  flow 
case  and  3.33,  3.34  for  the  unsteady  flow  case.  The  two  linear  equations^  are  then 
solved  again  with  Gaussian  elimination  to  obtain  the  vorticity  distributions.  For  the 
unsteady  case,  we  add  in  the  additional  requirement  of  finding  the  product  of  the 
tangential  velocities  at  the  upper  and  lower  trailing  edges.  These  are  demanded  to  be 
negative  as  there  are  basically  two  possible  solutions  to  the  Kutta  equations  due  to  its 
original  quadratic  nature. 

10.  Subroutine  NACA45 

This  subroutine  is  intended  to  calculate  the  camber  and  thickness  distribution 
of  the  NACA  4-digit  and  the  NACA  5-digit  airfoils  of  type  230XX  which  share  common 
thickness  distributions  with  the  4-digit  airfoil  having  the  same  thickness  to  chord  ratio. 
This  subroutine  is  called  by  subroutine  BODY  which  is  in  turn  called  by  subroutine 
SETUP  and  in  turn  called  by  the  MAIN  program. 

11.  Subroutine  NEWPOS 

This  is  a  new  subroutine  introduced  as  a  result  of  setting  up  the  five  frames  of 
reference.  Most  of  the  coordinates  computation  is  done  with  respect  to  the  respective 
local  frame  of  reference.  Its  purpose  then  is  to  transform  all  coordinates  in  the  respec- 
tive local  frames  of  reference  to  the  global  frame  of  referenced.  This  is  necessary  for  the 
computation  of  the  influence  coefficients  as  well  as  in  the  convection  of  the  wake  core 
vortices.  It  facilitates  the  simple  requirements  of  defining  the  airfoil  once  only  with  re- 
spect to  its  local  frame  of  reference.  Orientation  and  displacement  of  the  airfoil  in  the 
two-dimensional  plane  is  henceforth  calculated  through  this  subroutine.  This  subroutine 


16  For  the  unsteady  case,  the  original  two  equations  are  non-linear  and  were  subsequently 
linearised  in  the  discussion  in  Chapter  3. 

1 7  It  has  a  secondary  function  to  obtain  the  slopes  for  the  airfoil  panels  and  the  wake  element 
panels. 
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is  called  from  the  MAIN  program  in  several  locations  immediately  after  the  local  coor- 
dinates are  computed. 

12.  Subroutine  PRESS 

This  subroutine  calculates  the  pressure  distribution  over  the  airfoil  panels  after 
the  iterative  solution  for  the  unsteady  flow  problem  has  successfully  met  the  convergence 

criterion.    It  first  computes  the  tangential  velocities  at  all  panel  control  points  using 

P+  1 

equation  3.40  (with  p  =  — - — ),  then  performs  the  disturbance  potential  evaluation  at 

the  current  time  step  according  to  equations  3.41  through  3.43.  Together  with  the  dis- 
turbance potential  data  obtained  from  the  previous  time  step,  it  calculates  the  pressure 
distribution  using  equation  3.37. 

13.  Subroutine  SETUP 

This  subroutine  sets  up  the  local  panel  nodal  coordinates  for  MAIN  program 
by  reading  the  fourth  through  seventh  data  sets  of  the  input  file  if  IFLAG  =  1  is  set. 
It  skips  the  data  reading  if  IFLAG  =  0  and  proceeds  to  set  up  the  node  distribution  and 
calls  subroutine  BODY  to  calculate  the  airfoil  local  coordinates.  The  node  distribution 
adopts  a  cosine  formula  in  order  to  have  closely  packed  panels  toward  the  leading  and 
trailing  edges  for  improvements  in  solution  accuracy.  Regardless  of  how  the  nodal  co- 
ordinates are  obtained  ,  subroutine  SETUP  determines  the  local  panel  slopes  and  airfoil 
perimeter  length. 

14.  Subroutine  TEVVAK 

This  subroutine  is  intended  to  calculate  the  resultant  velocity  components  at  the 
mid-point  of  the  shed  vorticity  panel  with  respect  to  the  current  frozen  frame  of  refer- 
ence. These  velocity  components  are  necessary  to  ensure  that  the  correct  shed  vorticity 
panels'  length  and  orientation  are  established.  This  is  the  governing  criterion  for  the  it- 
erative solution  scheme  of  the  unsteady  flow  case.  This  subroutine  is  called  by  the 
MAIN  program  at  every  iteration  cycle  of  each  time  step  for  the  unsteady  flow  calcu- 
lation. 

15.  Subroutine  VELDIS 

This  subroutine  essentially  performs  the  same  function  as  subroutine  PRESS  for 
the  steady  flow  case.  It  is  thus  redundant  and  could  be  deleted  with  some  modification 
work  necessary  for  subroutine  PRESS.  This  work  has  been  implemented  by  Krainer  for 
the  single  airfoil  case.  This  subroutine  is  called  by  the  MAIN  program  in  the  steady  flow 
computation. 
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C.  INPUT  DATA  FOR  PROGRAM  USPOTF2 

The  Input  data  is  similar  to  that  for  the  single  airfoil  case.  Program  USPOTF2  reads 
its  input  data  from  filecode  1.  An  example  of  an  input  data  file  is  attached  in  Appendix 
B  for  the  case  when  the  airfoil  nodal  coordinates  are  input  by  the  user.  Computer  gen- 
erated airfoil  coordinates  are  another  option  that  can  be  selected  if  the  airfoil  chosen 
belongs  to  the  family  of  NACA  4-digit  or  5-digit  airfoils  of  type  230XX.  To  do  this,  the 
I  FLAG  parameter  is  set  to  zero  in  the  first  item  of  the  A'h  set  of  data  card  and  replace 
the  5th  and  &h  set  of  cards  by  two  cards  containing  the  particular  airfoil  NACA  number 
for  the  two  airfoils  using  format  (15).  Figure  7  contains  an  itemised  description  of  the 
sequential  input  variables. 

D.  OUTPUT  DATA  FROM  PROGRAM  USPOTF2 

The  Output  data  is  similar  to  that  for  the  single  airfoil  case.  Appendix  C  contains 
a  sample  output  data  generated  by  using  the  input  data  set  from  Appendix  B.  Due  to 
the  repetitive  nature  of  the  output  as  a  function  of  time,  only  selective  time  set  data  are 
shown.  The  output  data  file  begins  with  writing  out  what  the  program  has  read  from 
the  data  file  followed  by  the  computed  nodal  coordinates  only  if  they  are  generated  by 
the  program:  otherwise  it  proceeds  to  write  the  airfoil  computed  perimeter  length.  The 
next  set  of  output  data  are  the  steady  flow  solution  parameters  of  distributed  source 
strengths,  vorticity  strengths,  pressure-velocity  distribution,  force-moment  coefficient, 
potential  at  the  control  nodal  points  and  the  potential  at  the  leading  edge.  In  addition 
some  output  is  given  to  allow  an  assessment  of  the  accuracy  of  the  flow  solution.  This 
consists  of  the  normal  component  of  the  velocities  at  the  panel  control  points  and  the 
integral  of  the  disturbance  tangential  velocity  along  the  airfoil  contour.  If  only  the 
steady  flow  solution  is  required,  the  output  terminates  here. 

For  unsteady  flow,  in  addition  to  the  above,  the  unsteady  flow  parameters  are  also 
printed.  This  includes,  at  ev^ry  time  step,  the  iterative  printout  for  the  convergence  of 
the  length  and  orientation  of  the  wake  element  panel,  the  unsteady  solution  parameters 
on  the  airfoil  similar  to  that  of  the  steady  flow  parameters  as  well  as  the  trailing  wake 
core  vortices  data.  Again  the  same  checking  mechanisms  are  inserted  on  the  airfoil 
panels.  A  detailed  explanation  on  the  output  variable  names  is  listed  in  Figure  8.  All 
output  data  are  non-dimensional  quantities. 
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Data  Set  #1   Format  (15)  -  1  data  card. 

ITITLE        -  Number  of  title  cards  to  be  used  in  Data  Set  #2. 

Data  Set  #2   Format  (20A4)  -  ITITLE  data  cards. 
TITLE         -  Headings  printed  on  output  for  case  run 
identification. 


Data  Set  #3 

NAIRFO 

XSHIFT 

YSHIFT 


Data  Set  #4 
IFLAG 

NLOWER 
NUPPER 

Data  Set  #5 
x(I),y(I) 


Data  Set  #6 
x(I),y(I) 


Format  (I5,2F10. 6)  -  1  data  card. 

-  Number  of  airfoil,  in  this  case  =  2. 

-  Relative  X  distance  of  the  2  airfoil's  pivot  position 
with  respect  to  the  airfoil  global  coordinate  system. 

-  Relative  Y  distance  of  the  2  airfoil's  pivot  position 
with  respect  to  the  airfoil  global  coordinate  system. 

Format  (315)  -  1  data  card. 

-  0  if  airfoil  is  NACA  XXXX  or  230XX. 

-  1  otherwise. 

-  Number  of  panels  used  on  both  airfoil  lower  surface. 

-  Number  of  panels  used  on  both  airfoil  upper  surface. 

Format  (6F10. 6)  if  IFLAG  =  1  -  variable  data  cards. 

-  local  non-dimensionalised  x-nodal  followed  by  y-nodal 
coordinates  for  airfoil  1.   Total  of  Nlower+Nupper+1 
nodal  points  divided  into  6  points  per  data  card. 

Format  (6F10. 6)  -  variable  data  cards 

-  local  non-dimensionalised  x-nodal  followed  by  y-nodal 
coordinates  for  airfoil  2.   Total  of  Nlower+Nupper+1 
nodal  points  divided  into  6  points  per  data  card. 


Figure  7.   List  of  Input  Variables. 
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Data  Set  #7 

Format  (6F10. 6,13)  -  1  data  card 

ALPI(l) 

-  Initial  angle  of  attack  (AOA)  for  airfoil  1  in  degree. 

ALPI(2) 

-  Initial  angle  of  attack  (AOA)  for  airfoil  2  in  degree. 

DALP 

-  Absolute  change  in  AOA  in  degree  for  non  oscillatory 

motions. 

-  Maximum  amplitude  of  AOA  change  in  degree  for 

rotational  harmomic  motions. 

TCON 

-  Non-dimensional  rise  time  (V    tic)   of 

AOA  for  motion  involving  modif ied-ramp  change  in  AOA. 

FREQ 

-  Non-dimensional  oscillation  (coc/F^) 

for  harmonic  motions. 

PIVOT 

-  The  length  from  the  leading  edge  to  the  pivot  point 

for  the  local  system. 

IPHASE 

-  Flag  for  in-phase  and  out-of-phase  motion. 

-  0  out-of-phase  motion. 

-  1  in-phase  motion. 

Data  Set  #8 

Format  (8F10.6)  -  1  data  card. 

UGUST 

-  Magnitude  of  non-dimensional  gust  velocity  along 

global  x-direction. 

VGUST 

-  Magnitude  of  non-dimensional  gust  velocity  along 

global  y-direction. 

DELHX(l) 

-  Non-dimensional  trans lational  chordwise  amplitude 

for  airfoil  1 

DELHX(2) 

-  Non-dimensional  translational  chordwise  amplitude 

for  airfoil  2 

Figure  7  (cont'd) 
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Data  Set  #8 
DELHY(l) 

DELHY(2) 

PHASE(l) 

PHASE(2) 


Data  Set  #9 
TF 

DTS 


TOL 

TADJ 
SCLA 

NGIES 


(Cont'd) 

-  Non-dimensional  translational  transverse  amplitude 
for  airfoil  1 

-  Non-dimensional  translational  transverse  amplitude 
for  airfoil  2 

-  Phase  angle  in  degree  between  the  chordwise  and 
transverse  translational  oscillation  with  the  latter 
as  reference  for  the  first  airfoil. 

-  Phase  angle  in  degree  between  the  chordwise  and 
transverse  translational  oscillation  with  the  latter 
as  reference  for  the  second  airfoil. 

Format  (5F10. 6,15)  -  1  data  card. 

-  Final  non-dimensional  time  to  terminate  unsteady  flow 
solution. 

-  Starting  time  step  for  non-osc.  motions  if  TADJ  =  0.  0 

-  No.  of  computational  steps  per  cycle  for 
harmonic  motion. 

-  Baseline  time  step  size  for  all  motions  if  TADJ  ^  0. 0 

-  Tolerance  criterion  for  checking  the  convergence 
between  successive  iterations  of  {Uw)^   and  (Vj^)^. 

-  Factor  by  which  DTS  will  be  adjusted. 

-  Steady  lift  coefficient  for  the  single  airfoil  at  the 
specified  AOA 

-  Option  for  changing  the  unsteady  Kutta  condition  to 
satisfy  the  tangential  velocity  as  per  Geising's  case 

-  0  equal  pressure  at  the  trailing  edge  panels. 

-  1  equal  tgt  velocities  at  the  trailing  edge  panels. 


Figure  7  (cont'd) 
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TK       -  Time  step  t^. 

TKM1     -  Time  step  t^_^. 

ALPHA(L)  -  Angle  of  attack  of  airfoil  L  at  time  t^. 

OMEGA(L)  -  Rotational  velocity  (positive  counter  clockwise)  at  time 

tfc      for  airfoil  L 
U(L)     -  Chordwise  translational  velocity  (positive  forward)  at 

time  tfc   for  airfoil  L 
V(L)     -  Transverse  translational  velocity  (positive  downward)  at 

time  tfc   for  airfoil  L 
NITR     -  Iteration  number. 

VXV(L)  -  Iterative  solution  of  (U^)^  of  airfoil  L 
VYV(L)  -  Iterative  solution  of  iV^)^  of  airfoil  L 
WAKE(L)   -  Iterative  solution  of  shed  vorticity  panel  length  A^ 

of  airfoil  L 
THETA(L)  -  Iterative  solution  of  shed  vorticity  panel  orientation  0^ 

of  airfoil  L 
GAMK(L)   -  Iterative  solution  of  the  strength  of  the  current 

vorticity  distribution  of  airfoil  L 
J        -  Panel  number. 

XI(J)     -  local  x-coordinate  of  the  midpoint  of  jt"   panel. 
YI(J)     -  local  y-coordinate  of  the  midpoint  of  jt"   panel. 
X(J)     -  global  X-coordinate  of  the  midpoint  of  jt"   panel. 
Y(J)     -  global  Y-coordinate  of  the  midpoint  of  jt"   panel. 
Q(J)     -  Strength  of  source  distribution  on  the  j*-"   panel. 
CP(J)    -  Pressure  coefficient  at  the  midpoint  of  the  j*"   panel. 
V(J)     -  Total  tangential  velocity  at  the  midpoint  of  the  jt"   panel 

with  respect  to  the  moving  local  system. 


Figure  8.      List  of  Output  Variables. 
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VN(J) 

Norraal  velocity  at  the  midpoint  of  the  jt^1   panel 

with  respect 

to  the  moving  local  system. 

PHIK(J)   - 

Potential  at 

the  mid-point  of  the  jth   panel  at 

the  current  time  step. 

PHI(J) 

Potential  at 

the  mid-point  of  the  J™   panel  at 

previous  time  step. 

INTGAMMA  - 

Integral  of  the  disturbance  velocity  around  the  airfoil. 

CD(L) 

Drag  coefficient  of  airfoil  L. 

CL(L) 

Lift  coefficient  of  airfoil  L. 

CM(L) 

Pitching  moment  coefficient  about  leading  edge 

of  airfoil  L. 

M 

Trailing  wake  core  vortex  number. 

X1I(M) 

X-coordinate 

of  the  center  of  the  m^"   core  vortex 

of  airfoil  1 

with  respect  to  local  system. 

Y1I(M) 

Y-coordinate 

of  the  center  of  the  m*-"   core  vortex 

of  airfoil  1 

with  respect  to  local  system. 

X1(M) 

X-coordinate 

of  the  center  of  the  m^"   core  vortex 

of  airfoil  1 

with  respect  to  global  system. 

Y1(M) 

Y-coordinate 

of  the  center  of  the  m^   core  vortex 

of  airfoil  1 

with  respect  to  global  system. 

X2I(M) 

X-coordinate 

of  the  center  of  the  m^b   core  vortex 

of  airfoil  2 

with  respect  to  local  system. 

Y2I(M) 

Y-coordinate 

of  the  center  of  the  m*-"   core  vortex 

of  airfoil  2 

with  respect  to  local  system. 

X2(M) 

X-coordinate 

of  the  center  of  the  at"   core  vortex 

of  airfoil  2 

with  respect  to  global  system. 

Y2(M) 

Y-coordinate 

of  the  center  of  the  m^"   core  vortex 

of  airfoil  2 

with  respect  to  global  system. 

CIRC(M,L)- 

Circulation  strength  of  the  wi^h   core  vortex  of 

airfoil  L 

Figure  8  (Cont'd) 
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V.     RESULTS  AND  DISCUSSION  OF  CASE-RUNS 

USPOTF2  is  primarily  written  as  a  follow-up  to  U2DIIF  for  the  single  airfoil.  All 
the  case-runs  with  the  exception  of  the  gust  case  will  be  presented.  The  approach  in  the 
gust  case  is  not  consistent  with  the  requirements  for  an  irrotational  flow  field  and  will 
not  be  treated  in  this  report.  The  step  change  in  angle  of  attack  (AOA)  will  be  compared 
with  Giesing's  for  the  same  airfoil  set  at  the  final  AOA,  undergoing  an  impulsive  start 
from  rest.  As  there  exist  no  comparison  data  for  the  other  sub-cases,  the  results  will  thus 
not  be  as  extensive  as  the  step  input.  They  are  documented  here  for  the  sole  purpose 
of  illustrating  the  capability  of  the  code. 

A.     STEP  CHANGE  IN  ANGLE  OF  ATTACK 

1.  Case  Run  Definition 

Consider  a  first  case  of  two  airfoils  initially  at  zero  AOA  to  the  free  stream  V^ 
which  undergo  an  out-of-phase  step  change  in  AOA  (a„„)  at  time  i0.  Consider  a  second 
case  of  two  airfoils  at  rest,  set  initially  at  an  AOA  of  astep  out-of-phase  with  one  another, 
given  an  impulsive  start  to  \\c.  The  above  two  cases  are  equivalent  within  the  thin 
airfoil  approximation.  Case  1  is  computed  by  USPOTF2  and  case  2  is  obtained  by 
Giesing's  computational  analysis. 

2.  Differences  between  LSPOTF2  and  Giesing's  code 

While  both  codes  use  an  extension  of  the  PANEL  method  to  solve  for  the  un- 
steady potential  flow  solution  for  two  airfoils,  perfect  correlation  is  not  possible  for  se- 
veral reasons.  The  reader  is  referred  to  References  [5]  and  [7]  for  a  detailed  description 
of  Giesing's  approach.    Some  basic  differences  are: 

•  USPOTF2  models  the  wake  vortex  sheet  by  a  wake  element  and  a  series  of  point 
vortices  shed  through  the  wake  element.  This  follows  the  approach  of  Basu  and 
Hancock  with  the  necessary  assumptions  on  the  wake  element  characteristics. 
Giesing  treated  the  wake  vortex  sheet  to  comprise  of  a  distributed  line  vortex 
which  is  convected  at  the  local  fluid  velocity  at  the  vortex  location  assuming  that 
the  small  portion  of  the  wake  being  shed  does  not  contribute  to  the  convection  of 
the  wake. 

•  The  time  step  increment  in  LSPOTF2  is  a  parameter  that  affects  the  overran  results 
of  the  unsteady  flow  in  that  by  having  small  time  step,  the  point  vortices  being  shed 
will  be  greater  in  number  but  weaker  in  strength,  while  for  bigger  time  step,  the 
point  vortices  become  fewer  but  have  effectively  stronger  strengths.  The  time  step 
increment  in  Giesing's  code  is  important  in  that  the  assumption  that  the  small 
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portion  of  the  wake  being  shed  does  not  contribute  to  the  convection  of  the  wake 
is  only  exactly  true  in  the  limit  when  the  time  step  becomes  zero. 

Giesing  uses  the  Adam's  formula  of  varying  degreel8  to  convect  the  wake  vortex 
sheet  while  USPOTF2  uses  only  a  predictor  algorithm. 

The  treatments  of  the  circulation  T(/)  are  different  for  both  codes.  USPOTF2  ex- 
tended the  influence  coefficient  concept  to  the  unsteady  flow  regime  and  solved  for 
the  circulation  using  the  Kutta  condition  of  equations  3.33  and  3.34.  All  calcu- 
lations were  done  in  the  moving  frame  in  the  presence  of  the  unsteady  wake  for- 
mations. Giesing  considered  the  circulation  to  be  a  combination  of  the 
quasi-steady  circulation  and  the  circulation  due  to  the  vortex  wakes.  From  Refer- 
ence [  5] 

no  =  rQ([)  +  rs(f)  (5.i) 

where  ro(l)  is  the  circulation  required  to  satisfy  the  Kutta  condition  on  body  (/)  as 
it  moves  through  the  fluid  when  it  is  assumed  that  the  body  does  not  shed  any 
vorticity  and  Tr(/)  is  that  circulation  required  to  satisfy  the  Kutta  condition  on  body 
(/)  as  it  travels  through  the  flow  field  generated  by  the  vortex  wakes  shed  by  the  two 
bodies  and  the  flow  field  generated  by  the  circulatory  flow  about  the  other  body. 

Giesing's  Kutta  condition  reqires  equal  tangential  velocities  at  the  upper  and  lower 
surface  panels  at  the  trailing  edges  while  the  Kutta  condition  of  USPOTF2  pre- 
scribes equal  pressures.  In  order  to  have  a  meaningful  comparison,  LSPOTF2  in- 
cludes an  option  for  equal  tangential  velocities. 

USPOTF2  requires  the  actual  computation  of  the  total  velocity  potential  from  a 
combination  of  an  onset  flow  potential,  disturbance  potential  and  an  assumption 
of  a  reference  potential.  Giesing  uses  the  technique  of  the  Douglas  Potential  Flow 
program  which  treats  the  potential  as  the  combination  of  the  quasi-steady  potential 
and  the  potential  due  to  the  vortex  wakes  where  these  two  terms  are  defined  as 
before.  The  total  velocity  potential,  in  his  case,  need  not  be  computed  but  only  the 
time  derivative  of  it  which  is  then  written  in  terms  of  previously  known  parameters. 

The  number  of  panels  and  its  distribution  on  the  airfoils  are  different  for  both 
codes.  This  will  be  accentuated  at  the  trailing  edges  and  will  lend  itself  to  differ- 
ences in  the  trailing  edge  panel  distribution. 

3.     Results  and  Discussions 

Figure  9  shows  Giesing's  results  for  the  Von  Mises  8.4  per-cent  thick,  symmet- 
rical airfoil  undergoing  an  impulsive  start.  This  can  be  compared  with  Figure  10  which 
shows  the  result  obtained  with  USPOTF2  for  the  same  time  step  with  NGIES  set  to 
onei9.  The  results  are  not  perfectly  correlated  but  the  quality  and  order  of  magnitude 
agreement  are  excellent.   Figure  11  gives  essentially  the  same  result  but  with  NGIES  set 


• 
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1 8  This  is  essentially  a  predictor-corrector  algorithm. 

19  This  Kutta  condition  results  in  equal  tangential  velocities  at  the  trailing  edge  panels  of  the 
airfoils. 
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to  zero-0.  Surprisingly,  the  results  obtained  for  both  types  of  Kutta  condition  turn  out 
to  be  quite  similar.  This  is  seen  especially  in  the  pressure  coefficient  plots  when  the  two 
plots  are  put  together  as  seen  in  Figure  12.  The  pressure  coefficient  agrees  over  70  per- 
cent of  chord  length  at  the  lower  surface  and  over  90  per-cent  of  chord  length  at  the 
upper  surface  with  the  greatest  discrepancy  at  the  trailing  edge.  Further  comparison  for 
smaller  angles  of  attack  would  be  of  interest  to  see  whether  this  similarity  is  generally 
true. 

Figure  13  compares  the  aerodynamic  characteristics  when  the  vertical  distance 
VSHIFT  is  varied.  Figure  14  gives  the  time  variation  for  a  larger  time  step  of  the 
normalised  lift,  moment  coefficient  and  drag  coefficient  for  the  particular  case  of 
YSHIFT  =  2.0. 


20  This  Kutta  condition  results  in  equal  pressure  coefficients  at  the  trailing  edge  panels  of  the 
airfoils. 
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(a)  Pressure  Distribution  and  Vortex  locations. 


Figure  9.  Giesing's  Calculated  Results:  Impulsive  Start  for  an  8.4  per-cent  thick 
Von  Mises  airfoil  set  at  a  =  0.8  radians  for  YSHIF.T  =  2.0  [reproduced 
with  permission  from  Ref.  5]. 
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(a)  Pressure  Distribution  at  tl'Jc  =  0.65 

Figure  10.  USPOTF2  Results  obtained  with  USPOTF2  for  equal  velocities  at  the 
trailing  edge.:  Step  change  in  AOA  for  an  8.4  per-cent  thick  Von 
Mises  airfoils  placed  at  2  chord  length  vertical  distance  with  initial 
AOA  =  0.0  radians  and  final  AOA  =  0.8  radians  pivoting  at  the 
leading  edges. 
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Figure  10  (Cont'd) 
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Figure  10  (Cont'd) 
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(a)  Pressure  Distribution  at  tV^jc  =  0.65 


Figure   11.      ISPOTF2  Results  obtained  with  USPOTF2  for  equal  pressures  at  the 

trailing  edge.:  Step  change  in  AOA  for  an  8.4  per-cent  thick  Von 
Mises  airfoils  placed  at  2  chord  length  vertical  distance  with  initial 
AOA  =  0.0  radians  and  final  AOA  =  0.8  radians  pivoting  at  the 
leading  edges. 
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Figure  11  (Cont'd) 


62 


q 
CO 

6 

CO 

b 

o 

_i 
o 

* 
d 

CM 

d 

o 

(c)  Tin 

| | ] | r r 

)                    0.2                   0.4                   0.6                   0.8 

NON-DIMENSIONALISED  TIME 

ie  History  of  the  Lift  Coefficients. 

Figure  11  (Cont'd) 
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Figure  12.      Pressure  Distributions  with  different  Kutta  condition. 
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(a)  Pressure  Distribution  at  tl'Jc  =  0.65 


Figure   13.      Pressure  Distribution  and  Lift  History  as  a  function  of  the  Vertical  dis- 


tance between  the  airfoils. 
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(b)  Time  History  of  the  Lift  Coefficients. 


Figure  13  (Cont'd) 
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(a)  Time  Historv  of  the  Lift  Coefficients. 


Figure  14.       Step  Change  in  Angle  Of  Attack:     Step  change  in  AOA  for  an  8.4 

per-cent  thick  Von  Mises  airfoils  placed  at  2  chord  length  vertical  dis- 
tance with  initial  AOA  =  0.0  radians  and  final  AOA  =  0.8  radians 
pivoting  at  the  leading  edges. 
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(b)  Time  History  of  the  Moment  Coefficients. 


Figure  14  (Cont'd) 
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(c)  Time  History  of  the  Drag  Coefficients. 


Figure  14  (Cont'd) 
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(d)  Time  History  of  the  Circulation 


Figure  14  (Cont'd) 

B.     OTHER  SUB-CASES 

The  reader  is  referred  to  the  results  and  disussions  of  Reference  [  1  ]  for  the  com- 
parison of  the  single  airfoil  case  with  existing  codes  in  the  phase  and  magnitude  re- 
lationship of  the  aerodynamic  coefficients  to  the  forcing  function  as  well  as  the  wake 
convection.  The  present  results  show  that  the  same  trend  is  maintained  as  in  the  single 
airfoil  case  with  no  significant  phase  shift  but  with  a  general  magnitude  change  due  to 
an  effective  ground  effect  caused  by  the  presence  of  the  second  airfoil. 
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1.  Modified  Ramp  Change  in  Angle  of  Attack 

As  for  the  single  airfoil,  the  modified  ramp  is  defined  mathematically  as  follows: 

0  t  <  0 

ot(r)    =    6a{3  -  2//t)/2/t2     0  <  /  <  t  (5.2) 

(5a  t  >  x 

where  6y.  is  the  magnitude  of  the  AOA  change  and  t  is  the  rise  time  for  the  AOA  to 
reach  its  final  value.  Figure  15  treats  the  case  of  the  2  airfoils  undergoing  an  out-of- 
phase  modified  ramp  change  in  the  angle  of  attack  of  0.1  radians  with  a  time  constant 
of  1.5.  The  plots  show  the  effects  on  the  aerodynamic  coefficients  due  to  the  presence 
of  the  second  airfoil.  The  corresponding  aerodynamic  coefficients  of  the  single  airfoil 
are  plotted  for  comparison  purposes. 

2.  Translational  Harmonic  Motion 

The  code  is  capable  of  computing  the  unsteady  flow  solution  for  any  general 
translational  motion  described  by  a  chordwise  and  a  transverse  component  bearing  a 
given  phase  relationship  with  the  restriction  that  the  2  airfoils  move  only  in-phase  or 
out-of-phase.    The  translational  harmonic  motion  is  described  by 

hy{t)  =  dhy  Sin  {cot) 

hx(t)  =  6hx  Sin  {cor  +  /.)  (5.3) 

where  a  is  the  oscillation  frequency,  /  is  the  phase  angle  between  the  chordwise  and 
transverse  oscillation  and  dhx  and  dhy  are  the  magnitudes  of  chordwise  and  transverse 
oscillations  respectively.  The  case-run  considered  in  this  section  relates  to  a  pure  heav- 
ing or  plunging  motion.  A  NACA-0015  airfoil  is  chosen  for  the  case-run.  The  airfoils 
are  set  at  zero  radian  angles  of  attack  and  subsequently  given  an  out-of-  phase  plunging 
oscillation  at  an  amplitude  of  0.018  chord  length  at  a  non-dimensionalised  frequency  of 
1.  Figure  16  shows  the  aerodynamic  coefficients  for  the  two  airfoils  and  compares  them 
with  the  single  airfoil  case  where  applicable.  Note  that  the  plots  for  the  trailing  wake 
uses  different  scales  for  the  x  and  y  axes. 

3.  Rotational  Harmonic  Motion 

The  treatment  of  the  harmonic  pitching  motion  is  similar  to  the  modified  ramp 
case.  As  for  the  other  sub-cases,  the  airfoils  are  restricted  to  in-phase  and  out-of  phase 
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motion.  The  case  of  the  out-of-phase  motion  will  be  treated  here.  The  harmonic 
pitching  oscillation  is  described  by: 

o.U)  =  da  Sin  (cot)  (5.4) 

where  da,  and  a>  are  the  amplitude  and  frequency  of  the  harmonic  oscillation  respec- 
tively. Figure  17  shows  the  results  of  the  8.4%  thick  Von  Mises  symmetric  airfoil  os- 
cillating at  an  amplitude  of  0.1  radian  at  a  reduced  frequency  ofojcjV^  =  20.0  about  the 
leading  edge.  Again  the  plots  are  given  together  with  the  single  airfoil  case  undergoing 
the  same  motion  for  comparison  purposes. 
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(a)  Time  History  of  the  Lift  Coefficients. 


Figure  15.  Modified  Ramp  Change  in  Angle  Of  Attack:  Ramp  change  in  AOA 
for  an  8.4  per-cent  thick  Von  Mises  airfoils  placed  at  1  chord  length 
vertical  distance  with  initial  AOA  =  0.0  radians  and  final  AOA  =  0.1 
radians,  rise  time  of  1.5  pivoting  at  the  mid  chord. 
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Figure  16.  Harmonic  Plunging  Motion:  Translational  harmonic  plunging  AOA 
for  XACA-0015  airfoils  placed  at  1  chord  length  vertical  distance  set 
at  AOA  =  0.0  radians  with  plunging  amplitude  of  0.018  chord  length 
at  a  reduced  frequency  of  1. 
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Figure  17.  Harmonic  Pitching  Motion:  Rotational  harmonic  pitching  AOA  for 
XACA-0015  airfoils  placed  at  1  chord  length  vertical  distance  set  at 
AOA  =  0.0  radians  with  pitching  amplitude  of  0.1  radian  at  a  reduced 
frequency  of  4  pivoting  about  the  leading  edges. 
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VI.     CONCLUSION 

A.  GENERAL  COMMENTS 

USPOTF2  has  been  developed  as  an  intermediate  stage  to  getting  a  solution  for  a 
full  cascade  undergoing  unsteady  motion.  In  itself,  it  can  be  used  to  simulate  unsteady 
wing-tail,  wing-canard,  wing-flap,  wing-aileron,  horizontal  tail-elevator,  vertical 
tail-rudder  and  a  whole  host  of  moving-stationary  airfoil  interactions.  These  simu- 
lations require  a  little  addition  to  the  main  program  and  would  appear  as  modules 
modelling  the  variation  of  the  global  angles  of  attack  and  relative  displacements  as  per 
all  the  sub-cases  done.  The  subroutines  would  not  be  affected  by  the  severity  of  the  test 
case  except  for  possibly  subroutine  NEWPOS. 

Validation  of  LSPOTF2  has  been  done  against  Giesing's  code  for  the  particular  case 
of  a  step  change  in  angle  of  attack.  This  was  shown  to  have  good  correlation.  However, 
it  can  not  be  said  that  the  agreement  will  hold  for  other  angles  of  attack  and  more  sub- 
case runs  for  smaller  angles  of  attack  should  be  made.  In  the  same  manner,  the  success 
of  the  step  change  in  AOA  in  no  way  validates  the  other  sub-cases  which,  for  the  time 
being,  remain  unproven. 

B.  ENHANCING  USPOTF2  PROGRAM  CAPABILITY 

As  noted  above,  USPOTF2  needs  to  be  tried  more  extensively  either  with  the  exist- 
ing sub-cases  or  with  new  sub-cases  against  existing  numerical  or  experimental  results. 

The  software  in  itself  has  an  implicit  weakness  in  the  numerical  computation  of  the 
velocity  potential.  The  velocity  potential  at  the  leading  edge  is  obtained  by  integrating 
the  velocity  field  from  the  leading  edge  to  a  point  100  chord  lengths  upstream  of  the 
leading  edge.  Two  assumptions  are  made:  First,  the  disturbance  velocity  at  the  point 
100  chord  lengths  of  the  leading  edge  approaches  zero.  Second,  the  contribution  to  the 
velocity  potential  due  to  the  integration  from  the  point  100  chord  lengths  upstream  of 
the  leading  edge  to  upstream  infinity  must  not  change  in  time.  The  coefficients  of  pres- 
sure will  be  accurate  only,  if  those  two  assumptions  hold  (at  least  in  an  approximate 
sense).  This  however  does  not  give  a  very  'exact'  solution  as  by  increasing  the  chord 
length  by  900  per-cent  to  1000  chord  lengths,  there  is  a  change  in  the  pressure  coefficient 
at  the  trailing  edge  of  about  8  per-cent  for  the  first  time  step.  Krainer  has  implemented 
an  analytical  solution  to  the  potential  problem  for  the  single  airfoil.  A  similar  procedure 
to  upgrade  LSPOTF2  should  be  considered. 
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Other  improvements  to  the  code  would  be  the  reduction  of  redundant  computations. 
An  obvious  example  would  be  the  subroutines  VELD1S  and  PRESS  which  essentially 
do  the  same  work  for  the  steady  and  unsteady  flow  respectively.  Again  Krainer  [Ref. 
6]  has  improved  the  original  code  for  the  single  airfoil  (U2DIIF)  and  though  some  of 
his  improvements  were  implemented  in  USPOTF2,  there  still  remains  a  task  to  do  the 
full  job  completely. 

Finally,  the  original  primary  objective  of  extending  the  code  to  solve  for  the  un- 
steady flow  solution  for  3,4  airfoils  and  leading  to  a  full  cascade  still  remains  a  big  task, 
not  just  for  the  programmer  but  also  for  the  computer  system  in  terms  of  computational 
storage  and  time  requirements. 
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APPENDIX      A.  USPOTF2  SOURCE  LISTINGS 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


PROGRAM     USP0TF2 

VERSION      1 

SEPTEMBER  88 

UNSTEADY  MOTION  FOR  TWO  AIRFOILS  IN  POTENTIAL 

INCOMPRESSIBLE  FLOW 
USING  PANEL  METHODS  BASED  ON  THE  HESS  &  SMITH 
AND  WAKE  MODEL  BASED  ON  BASU  &  HANCOCK 


C 
C 
C 
C 
C 
C 
C 

c 
c 
c 


+ 
+ 
+ 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y(202) , 

COSTHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 
NP5,XSHIFT,YSHIFT,NAIRFO,XI(202),YI(202), 
COSTHL(201),SINTHL(201) 
/WAK/  VYW(2),VXW(2),WAKE(2),DT 
/WAK2/  VYWK(2) ,VXWK(2) 
/SING/  Q(200),GAMMA(2),QK(200),GAMK(2) 
/CORV/  CV(2,200),XC(2,200),YC(2,200),M,TD,CVVX(2,200), 

CVVY(2,200),XCI(2,200),YCI(2,200) 
/POT/  PHI(200),PHIK(200) 
/GUST/  UG(200),VG(200),XGF,UGUST,VGUST 
/EXTV/  UE(200),VN(200) 

/PARD/  NACAD(2),TAUD(2),EPSMAD(2),PTMAXD(2) 
/COF/  A(201,211),KEQNS 
/GEOM/  SINALF(2),COSALF(2),OMEGA(2),UX(2),UY(2), PIVOT, 

XPRM,YPRM 
/CPD/  CP(200),SCL,T,SCM,SGAM 
/GIES/  NGIES 
DIMENSION  XXC(2,200),YYC(2,200),TOL1(2),TOL2(2),THENP1(2), 
+SINANG(2),COSANG(2),COSDA(2),SINDA(2),DHX(2),DHY(2),ALP(2), 
+ALPHA(2) ,ALPI(2) ,ANGLE(2) ,DELHX(2) ,DELHY(2) ,PHASE(2) ,PHA(2) , 
+HX(2),HXO(2),HY(2),HYO(2) 


COMMON 
COMMON 
COMMON 
COMMON 
h 

COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
h 

COMMON 
COMMON 


c 

c 

INITIALISATION 

c 

PI      =  3. 1415926585 

M       =0 

XPRM    =  0.  0 

YPRM    =  0.  0 

c 

c 

INPUT  FROM  FILE  CODE  5 

c 

AND  SET  UP  LOCAL  PANEL  NODES  AND  SLOPES 

WRITE  (6,1003) 
1003  FORMAT  (/////,'  DATA  READ  FROM  FILE  CODE  l',//) 

CALL  INDATA 

CALL  SETUP 

READ  (1,502)  ALPI(l),ALPI(2),DALP,TCON,FREQ,PIVOT,IPHASE 

WRITE  (6,502)  ALPI( 1) ,ALPI(2) ,DALP,TCON,FREQ,PIVOT, IPHASE 
501   FORMAT  (8F10. 6) 
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502  FORMAT  (6F10. 6,13) 

READ  (1,501)  UGUST,VGUST,DELHX(1),DELHX(2),DELHY(1),DELHY(2), 
+  PHASE(1),PHASE(2) 

WRITE  (6,501)  UGUST,VGUST,DELHX(1),DELHX(2),DELHY(1),DELHY(2), 
+  PHASE(1),PHASE(2) 

READ  (1,503)  TF,DTS,TOL,TADJ,SCL,SCM,SGAM,NGIES 
WRITE  (6,503)  TF,DTS,TOL,TADJ,SCL,SCM,SGAM,NGIES 

503  FORMAT  (7F10.6,I3) 

IF  (IFLAG  . EQ.  0)  WRITE  (6,1005) 
1005  FORMAT  (///,'  COORDINATES  OF  AIRFOIL  NODES', 
+  //,3X,'  X/C,6X,'  Y/C',/) 
IF  (IFLAG  .EQ.  0)  WRITE  (6,1010)  (XI( I) ,YI( I) , I=l,2*(NODTOT+l)) 
1010  FORMAT  (F10.  6.F10.  6) 

WRITE  (6,1000)  NAIRFO 
1000  FORMAT(//,' TOTAL  NO  OF  AIRFOILS  =  ',14,/) 
C 


c 

STEADY 

FLOW  CALCULATION  A 

c 

c 
c 

INITIAL 

DEFINITION 

NP1 

=  NODTOT  +  1 

NP2 

=  2  *  NP1  +  1 

NP3 

=  NAIRFO*NODTOT+l 

NP4 

=  NP3+1 

NP5 

=  NP4+1 

DO  101  I  =  1,2 

ALPHA(I)  =  ALPI(I) 

IF  (ALPHA(I)  .GT.  90.)    GO  TO  200 

COSALF(I)  =  C0S(ALPHA(I)*PI/180. ) 

101  SINALF(I)  =  SIN(ALPHA(I)*PI/180.  ) 
CALL  NEWPOS(O) 

DO  1100  L  =  1, NAIRFO 
1100  WRITE  (6,1020)  L,SS(L) 

WRITE  (6,1040)  XSHIFT,YSHIFT 
1020  FORMAT(//,'  AIRFOIL(  '  12 , '  )  PERIMETER  LENGTH  =  \F10.6,/) 
1040  FORMAT(//,'  XSHIFT  =  f,F5.1,'    YSHIFT  =  \F5.1,/) 

DO  102  I  =  1,2 

102  WRITE  (6,1030)  I,ALPHA(I) 

1030  FORMAT  (//,'  STEADY  FLOW  SOLUTION  AT  ALPHA( ' , 12 , ' )  =  '  ,F10.6,/) 
IF  (SCL. NE.  0)   WRITE  (8,1035) 

1035  FORMAT(3X,'TIME' ,4X, ' CL( 1)/SCL' , IX, ' CL(2)/SCL' , IX, *CM( 1)/SCM' , 
+1X,'CM(2)/SCM' ,3X,'CD(1)' ,4X,'CD(2)' ,4X, 

+'GAMK(1)/SG' ,1X,*GAMK(2)/SG' ,//) 
IF  (SCL. EQ. 0)   WRITE  (8,1036) 

1036  FORMAT(3X,'TIME' ,6X,'CL(1)' ,5X, 'CL(2) ' ,5X, 'CM( 1) ' , 
+5X,'CM(2)' ,5X,'CD(1)' ,4X, ' CD(2) ' //) 

CALL  INFL  (0) 
CALL  COEF  (0) 
CALL  GAUSS(3,0,0) 
CALL  KUTTA(NITR,PVTAG) 
CALL  VELDIS 
SIN1  =  SINALF(l) 
SIN2  =  SINALF(2) 
COS1  =  COSALF(l) 
COS2  =  COSALF(2) 
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CALL  FANDM(SIN1,SIN2,C0S1,C0S2) 


c 
c 
c 

INITIALISATION  FOR  UNSTEADY  FLOW  CALCULATION 

DA 

0.0 

XGF 

0.0 

DO  100  L 

=  l,NAIRFO 

HX(L) 

0.  0 

HY(L)   = 

0.  0 

HXO(L)   = 

0.0 

HYO(L)   = 

0.0 

ANGLE (L) 

=  ALPI(L)*PI/180.  +  ATAN(VGUST/( 

OMEGA(L) 

=  0.  0 

ALP(L) 

=  ALPI(L) 

DHX(L) 

=  0.0 

DHY(L) 

=  0.  0 

COSANG(L) 

=  COS ( ANGLE (L)) 

SINANG(L) 

=  SIN(ANGLE(L)) 

UX(L) 

=  0.0 

UY(L) 

=  0.  0 

COSDA(L) 

=  1.  0 

SINDA(L) 

=  0.0 

KIG 

=  (L  -1)*N0DT0T 

VXW(L)  =  ( 

:OSALF(L) 

VYW(L)  =  ! 

SINALF(L) 

GAMK(L) 

=  GAMMA(L) 

PHA(L)   = 

PHASE(L)*PI/180. 

DO  100   IG  =  l,NODTOT 

UG(IG+KIG)   =0.0 

100 

VG(IG+KIG)   =  0.  0 

T 

0.0 

c 
c 
c 

TOLD 

0.0 

RIGID  ] 

30DY  MOTIONS  OF  AIRFOIL 

IF  (FREQ 

NE.  0.0)  GO  TO  1 

IF  (DALP 

EQ.  0.  0)  GO  TO  2 

IF  (TCON 

NE.  0.  0)  GO  TO  3 

IF  (IPHASE  .EQ.  0)  GO  TO  4 

ALPHA(1)= 

ALPI(l)  +  DALP 

ALPHA(2)= 

ALPI(2)  +  DALP 

GO  TO  5 

4 

ALPHA(1)= 

ALPI(l)  +  DALP 

ALPHA(2)= 

ALPI(2)  -  DALP 

5 

DO  6  I  = 

1,2 

COSALF(I) 

=  C0S(ALPHA(I)*PI/180. ) 

SINALF(I) 

=  SIN(ALPHA(I)*PI/180. ) 

6 

CONTINUE 

CALL  NEWPOS(O) 

3 

DT 

DTS 

TD 

DTS 

GO  TO  60 

2 

IF  ((UGUST  .EQ.  0.0)  .AND.  (VGUST  .  EQ.  0.0)) 

DT 

DTS 

TD 

DTS 

GO  TO  60 

TO  BEGIN 


GO  TO  200 
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1 

60 
105 

40 


DT  =  2. 0*PI/(FREQ*DTS) 

TD  =  DT 

T  =  DT 

WRITE  (6,1051) 

1  FORMAT  (/////,  '  ,*k**1Htic****1rk*1c***1r*****1tic4t4rk****ic****-it '  / 

+  '  ***  BEGIN 'UNSTEADY  FLOW  SOLUTION  ****«>/j 

+  '  *************************************** '  ) 

TRY     =0 

M       =  M  +  1 

IF  (T  .GT.  TF)  GO  TO  200 

STORE  CORE  VORTEX  COORDINATES  FOR  TIME  STEP  ADJUSTMENTS 


51 

50 


33 


41 


42 


43 


IF  (M  .EQ. 
DO  51   L 


1)  GO  TO  50 
=  l,NAIRFO 


DO  51    I  =  l.M-1 
XXC(L,I)   =  XCI(L,I) 
YYC(L,I)   =  YCI(L,I) 
IF  (FREQ  .NE.  0.  0)  GO  TO  11 
IF  (DALP  .EQ.  0.0)  GO  TO  22 
IF  (TCON  .NE.  0.  0)  GO  TO  33 

STEP  CHANGE  IN  AOA 

IF  (TADJ  .NE.  0.0)  GO  TO  70 

IF  INCREMENTAL  PROGRESSIVE  TIME  STEP  IS  REQUIRED 

USE  ...  TD  =  FLOAT(M+l)*DTS  ...  OTHERWISE  CONSTANT  TIME  STEP 

TD      =  DTS 

GO  TO  70 

MODIFIED  RAMP  CHANGE  IN  AOA 

IF  (T  .GT.  TCON)  GO  TO  34 

DAL     =  DALP  *  (3. -2. *T/TCON)*(T/TCON)**2 

OMEGA(l)    =  -  (DALP---P 1/180.  )  *  ( 6.  *T/(TCON*TCON)  )  *  (l.-T/TCON) 
IF  (IPHASE  .EQ.  0)  GO  TO  41 
ALPHA(1)=  ALPI(l)  +  DAL 
ALPK2)  +  DAL 
=  OMEGA(l) 


ALPHA(2)= 

0MEGA(2) 

GO  TO  42 

ALPHA(1)= 

ALPHA(2)= 

0MEGA(2) 

DO  43   I  = 

COSALF(I) 

SINALF(I) 

DA 

COSDA(I) 

SINDA(I) 

DHX( I ) 

DHY(I) 

UY(I) 

CONTINUE 

MTCON   =  M 

CALL  NEWPOS(O) 

GO  TO  70 


ALPI(l)  +  DAL 

ALPI(2)  -  DAL 
=  -  OMEGA(l) 
1,2 
=  COS(ALPHA(I)*PI/180. ) 
=  SIN(ALPHA(I)*PI/180. ) 
=  ALPHA(I)  -  ALP(I) 
=  COS(DA*PI/180. ) 
=  SIN(DA*PI/180. ) 
=  PIVOT  *  (1. -COSDA(I)) 
=  -  PIVOT  *  SINDA(I) 
=  PIVOT  *  OMEGA(I) 
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34 

DAL 

0.0 

IF  (IPHASE  . EQ.  0)  GO  TO  45 

ALPHA(1)= 

ALPI(l)  +  DALP 

ALPHA(2)= 

ALPI(2)  +  DALP 

GO  TO  46 

45 

ALPHA(l) 

=  ALPI(l)  +  DALP 

ALPHA(2) 

=  ALPI(2)  -  DALP 

46 

DO  44  I  = 

1,2 

COSALF(I) 

=  COS(ALPHA(I)*PI/180.  ) 

SINALF(I) 

=  SIN(ALPHA(I)*PI/180.  ) 

DA 

=  0.  0 

COSDA(I) 

=  1.0 

SINDA(I) 

=  0.0 

OMEGA(I) 

=  0.0 

DHX( I ) 

=  0.0 

DHY(I) 

=  0.0 

UX(I) 

=  0.0 

UY(I) 

=  0.0 

44 

CONTINUE 

CALL  NEWPOS(O) 

IF  (TAD J  . 

,  NE.  0.  0)  GO  TO  70 

TD 

FL0AT(M+1-MTC0N)*DTS 

C 
C 

GO  TO  70 

SHARP  ] 

EDGE  GUST  (UGUST  and/or 

C 
C 

INCONSISTENT  ASSUMPTIONS  REQ 

C 
22 

XGF 

T 

DO  113  L 

=  1,NAIRF0 

LIG 

(L-1)*N0DT0T 

KIG 

(L-l)-'-NPl 

DO  110   IG  =  l,NODTOT 

111 


120 


121 


110 


VGUST)  NOT  PROVEN  DUE  TO 
UIRING  ROTATIONAL  FLOW 


UG(IG+LIG)   =  0.  0 
VG(IG+LIG)   =0.0 
XG      =  X(IG+KIG) 
XGP1    =  X(IG+KIG+1) 
IF  (IG  .LT.  NLOWER+1)  GO  TO 
IF  (XGF  .LE.  XG)  GO  TO  110 
IF  (XGF  .GE.  XGP1)  GO  TO  111 
FAC     =  (XGF  -  XG)/(XGP1  - 


120 


XG) 


UG(IG+LIG) 
VG(IG+LIG) 
GO  TO  110 
UG(IG+LIG) 
VG(IG+LIG) 
GO  TO  110 
IF  (XGF  .LE 
IF  (XGF  .GE 


=  UGUST*FAC 
=  VGUST*FAC 

=  UGUST 
=  VGUST 


XGP1)  GO  TO  110 
XG)  GO  TO  121 

FAC     =  (XGF  -  XGP1)/(XG  -  XGP1) 

UG(IG+LIG)   =  UGUST*FAC 

VG(IG+LIG)   =  VGUST*FAC 

GO  TO  110 

UG(IG+LIG)   =  UGUST 

VG(IG+LIG)   =  VGUST 

CONTINUE 


96 


113 


11 


131 


12 


141 


142 


143 


70 


CONTINUE 

IF  (XGF  .LE.  COSALF(L))  MGUST  =  M 

IF  (TADJ  .NE.  0.0)  GO  TO  70 

IF  (XGF  .GT.  COSALF(L))  TD  =  FLOAT(M+l-MGUST)*DTS 

GO  TO  70 

TRANSLATION  HARMONIC  OSCILLATION 


(DALP  .NE.  0.  0)  GO  TO  12 
I  =  l.NAIRFO 
=  DELHX(I)  *  SIN(FREQ*T  +  PHA(I)) 
=  DELHY(I)  *  SIN(FREQ*T) 
=  HX(I)  -  HXO(I) 
=  HY(I)  -  HYO(I) 

=  DELHX(I)*FREQ*COS(FREQ*T+PHA(I)) 
DELHY(I)*FREQ*COS(FREQ*T) 


IF 

DO  131 

HX(I)   = 

HY(I)   = 

DHX(I)   = 

DHY(I)   = 

UX(I)    = 

UY(I)    = 

CONTINUE 

XPRM    =  DHX(2) 

YPRM    =  DHY(2) 

CALL  NEWPOS(O) 

GO  TO  70 


DHX(l) 
DHY(l) 


ROTATIONAL  HARMONIC  OSCILLATION 

DAL     =  DALP*SIN(FREQ*T) 

OMEGA(l)   =  -  (DALP*PI/180. )  *  FREQ  *  COS(FREQ'T) 
IF  (IPHASE  . EQ.  0)  GO  TO  141 
ALPHA(1)=  ALPI(l)  +  DAL 
ALPI(2)  +  DAL 
=  OMEGA(l) 


ALPHA(2)= 

OMEGA(2) 

GO  TO  142 

ALPHA(1)= 

ALPHA(2)= 

OMEGA(2) 

DO  143   I 

COSALF(I) 

SINALF(I) 

DA 

COSDA(I) 

SINDA(I) 

DHX(I) 

DHY(I) 

UY(I) 

CONTINUE 

CALL  NEWPOS(O) 


ALPI(l)  +  DAL 

ALPI(2)  -  DAL 
=  -  OMEGA(l) 

=  1,2 

=  C0S(ALPHA(I)*PI/180. ) 
=  SIN(ALPHA(I)*PI/180.  ) 
=  ALPHA(I)  -  ALP(I) 
=  COS(DA*PI/180. ) 
=  SIN(DA*PI/180. ) 
=  PIVOT  *  (1.  -COSDA(I)) 
=  -  PIVOT  *  SINDA(I) 
=  PIVOT  *  OMEGA(I) 


TRANSFORM  CORE  VORTEX  COORDINATES  W.  R.  T.  NEW  AIRFOIL  POSITION 


IF  (M  .EQ. 

DO  85   L 

DO  90   I 

XCI(L,I) 

YCI(L.I) 

XCO 

YCO 

XCI(L,I) 


1)  GO  TO  80 
=  l,NAIRFO 
=  1,M-1 

=  XXC(L.I)  + 

=  YYC(L,I)  + 
XCI(L,I) 
YCI(L.I) 

=  XCO*COSDA(L)  - 


CWX(L,I) 
CWY(L,I) 


DT 
DT 


YCO*SINDA(L)  +  DHX(L) 
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90   YCI(L,I) 
85   CONTINUE 

CALL  NEWP0S(3) 
80   CONTINUE 

WRITE  (6,1001) 
1001  FORMAT  (/////, 

WRITE  (6,1004) 


=  XCO*SINDA(L)  +  YCO*COSDA(L)  +  DHY(L) 


1004  FORMAT  (/,'  ALPHA(l) 
+  '  OMEGA(l)   =  '  ,F10.6 

+/,'  U(D      =  ',F10 

+F10.6,/,'  V(l) 

+  IX,'  NITR   VXW(l) 


T,DT 

TIME  STEP  TK  =  ' ,F10. 6 , 10X, fTK  -  TKM1  =  ' ,F10. 6,/) 
ALPHA( 1) ,ALPHA(2) ,OMEGA( 1) ,OMEGA(2) ,UX( 1) ,UX(2) , 
UY(1),UY(2) 

",F10.6,/, 
5X,  " 


,F10. 6,5X,'  ALPHA(2)   = 
f  OMEGA(2)   =  ' ,F10. 6, 

6,5X,'   U(2) 

',F10.6,5X,'   V(2)     =  *  , 


VYW(l) 


+  VXW(2) 


VYW(2)   WAKE(2) 


WAKE(l) 
THETA(2) 


THETA(l) 
GAMK(2) 


F10.6,///, 

GAMK(l) 
',/) 


CALCULATE  THE  TRAILING  EDGE  WAKE  ELEMENT 


10 


NUM 
NITR 
DO  15  L 
KI 

WAKE(L)  = 
THENPl(L) 


=  l,NAIRFO 
=  (L-l)*(NODTOT+l) 
SQRT(  VYW(  L)*VYW(  L)+VXW(  L)*VXW(  L)  )*DT 
=  ATAN2(VYW(L),VXW(L)) 
C0STHL(NP1+KI)  =  C0S(THENP1(L)) 
15   SINTHL(NP1+KI)  =  SIN(THENP1( L) ) 

WRITE  (6,1002)  NITR,VXW(1),VYW(1),WAKE(1),THENP1(1),GAMK(1) 
+VXW(2),VYW(2),WAKE(2),THENP1(2),GAMK(2) 
1002  FORMAT  (  15  ,4F10.  6  ,E14.  6  ,4F10.  6  ,E14.  6) 

XI(NP2)   =  XI(NPl)  +  WAKE(1)*C0STHL(NP1) 
YI(NP2)   =  YI(NPl)  +  WAKE(1)*SINTHL(NP1) 
XKNP2+1)   =  XI(NPl)  +  WAKE(2)*C0STHL(2*NP1) 
YI(NP2+1)   =  YI(NPl)  +  WAKE(2)*SINTHL(2*NP1) 
CALL  NEWPOS(l) 
CALL  INFL  (NITR) 
CALL  COEF  (NITR) 
CALL  GAUSS ( 3, M, NITR) 
C     WRITE  (6, *)  *A(I,NP)' , ( A( I , 101) , 1=1 , 100) 
C     WRITE  (6,*)  fA(I,NTOT)f , (A( I , 102) , 1=1, 100) 
C     WRITE  (6,*)' BEFORE  KUTTA' 
CALL  KUTTA(NITR,PVTAG) 
IF  (PVTAG  .LT.  0.01)  GOTO  13 
NUM  =  NUM  +  1 
C     WRITE  (6,*)' NUM  =' ,NUM 
IF  (NUM  .GT.  1)  GOTO  13 
DO  7  I  =  1,2 


13 


24 


VXW(I)  = 

=  1 

VYW(I)  = 

=  0 

GOTO  10 

CALL  TEWAK 

DO  24  L 

=  1,NAIRF0 

TOLl(L) 

=  ABS(VYW(L) 

TOL2(L) 

=  ABS(VXW(L) 

VYW(L) 

=  VYWK(L) 

VXW(L) 

=  VXWK(L) 

-  VYWK(L))/VYWK(L) 

-  VXWK(L))/VXWK(L) 


IF  ((TOLl(l)  .LT.  TOL)  .AND.  (TOL2(l)  .  LT.  TOL) 
+.AND.  (T0L1(2)  .LT.  TOL)  .AND.  (TOL2(2)  .  LT.  TOL))  GOTO  20 


IF  (NITR  .  GT.  25)  THEN 

TOL     =  TOL*10 

WRITE  (6,1023)  TOL 

END  IF 

IF  (NITR  .GT.  50)  STOP 

NITR    =  NITR  +  1 

GO  TO  10 
20   WRITE  (6,1011)  NITR 

1011  FORMAT  (//,'  CONVERGED  SOLUTION  OBTAINED  AFTER  NITR  =  ',13) 
1023  FORMAT  (//,'  *****  TOLERANCE  CRITERIA  CHANGED  :  TOL  =  * ,F10. 6) 

CALL  PRESS 

IF  ((UGUST  .EQ.  0.0)  .AND.  (VGUST  .  EQ.  0.0))  GO  TO  300 

SIN1  =  SINANG(l) 

SIN2  =  SINANG(2) 

COS1  =  COSANG(l) 

COS2  =  COSANG(2) 

CALL  FANDM(SIN1,SIN2,C0S1,C0S2) 

GO  TO  400 
300   SIN1  =  SINALF(l) 

SIN2  =  SINALF(2) 

COS1  =  COSALF(l) 

COS2  =  COSALF(2) 

CALL  FANDM(SIN1,SIN2,C0S1,C0S2) 
400   CONTINUE 
C 

C     ADJUST  TIME  STEP  (TADJ  .  NE.  0.0)  IF  NECESSARY 
C 

IF  (TADJ  .EQ.  0.  0)  GO  TO  95 

WRITE  (5,2001) 
2001  FORMAT  (//,'  DO  YOU  WANT  TO  ADJUST  TIME  STEP  ?  0  -  NO,  1  -  YES') 

READ  (5,*)  IDT 

IF  (IDT  .EQ.  0)  GO  TO  95 

DT      =  TADJ  *  DT 

T       =  TOLD  +  DT 

WRITE  (6,1006) 
1006  FORMAT  (//,'  BACK-TRACK  COMPUTATION  AND  ADJUST  TIME-STEP',//) 
C 

C        WAKE  ELEMENT  LEAVES  TRAILING  EDGE  AS  A  CORE-VORTEX 
C 

95  DO  96   L  =  l,NAIRFO 

CV(L,M)    =  SS(L)*(GAMMA(L)-GAMK(L)) 
CVVX(L,M)   =  VXW(L) 

96  CVVY(L,M)   =  VYW(L) 

XCI(l.M)    =  XI(NPl)  +  0.  5*WAKE(1)*C0STHL(NP1) 
YCI(1,M)    =  YI(NPl)  +  0. 5*WAKE(1)*SINTHL(NP1) 
XCI(2,M)    =  XI(NPl)  +  0.5*WAKE(2)*COSTHL(NP1*2) 
YCI(2,M)    =  YI(NPl)  +  0.5*WAKE(2)*SINTHL(NP1*2) 
CALL  NEWPOS(2) 
WRITE  (6,1052) 
1052  FORMAT  (//,'  TRAILING  VORTICES  DATA',//, 

+  4X,'M'  ,4X,'X1I(M)'  ,5X,'Y1I(M)'  ,6X,'X1(M)'  ,6X,*Y1(M)'  , 

+  5X/CIRC11  ,6X,'X2I(M)'  ,5X,'Y2I(M)'  , 

+  5X,'X2(M)' ,6X,'Y2(M)' ,6X,'CIRC2' ,/) 
DO  900      I   =   1,M 
900     WRITE   (6,1050)    I ,XCI( 1 , I) ,YCI( 1 ,1) ,XC( 1,1) ,YC( 1,1) ,CV( 1 , I) , 

+  XCI(2,I),YCI(2,I),XC(2,I),YC(2,I),CV(2,I) 
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1050  FORMAT(I5,10F11.6) 

CALL  CORVOR 
C 

C        RE-INITIALISE  PARAMETERS  FOR  NEXT  TIME  STEP  CALCULATION 
C 

DO  30   L  =  l,NAIRFO 

LI      =  (L-l)*NODTOT 

HXO(L)   =  HX(L) 

HYO(L)   =  HY(L) 

GAMMA(L)=  GAMK(L) 

ALP(L)   =  ALPHA(L) 

DO  30   I  =  l.NODTOT 

PHKI+LI)   =  PHIK(I+LI) 
30   CONTINUE 

TOLD    =  T 

DT      =  TD 

T       =  T  +  TD 

GO  TO  40 
200  WRITE  (6,1024)  TOL 
1024  FORMAT  (IX, /•*****  TOLERANCE  CRITERION  USED:  TOL  * ,F10. 6) 

STOP 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C     SUBROUTINE  INDATA  C 

C  C 

C  SET  PARAMETERS  OF  BODY  SHAPE                      C 

C  FLOW  SITUATION,  AND  NODE  DISTRIBUTION              C 

C  C 

C  USER  MUST  INPUT                                   C 

C  NLOWER  =  NUMBER  OF  NODES  ON  LOWER  SURFACE      C 

C  NUPPER  =  NUMBER  OF  NODES  ON  UPPER  SURFACE      C 

C  PLUS  DATA  ON  BODY  AND  SUBROUTINE  BODY              C 

C  C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C     NACAD(I)  NACA  NUMBERS  FOR  THE  TWO  AIRFOILS 

C     TAUD(I)   MAX  THICKNESS  FOR  THE  TWO  AIRFOILS 

C     EPSMAD(I)  ...  MAX  CAMBER  FOR  THE  TWO  AIRFOILS 

C     PTMAXD(I)  ...  CHORDWISE  POSN  FOR  MAX  CAMBER  FOR  THE  TOO  AIRFOILS 

C     

SUBROUTINE  INDATA 
DIMENSION  TITLE(20) 

COMMON  /BOD/  IFLAG, NLOWER, NUPPER, NODTOT,X( 202) ,Y(202) , 
+  COSTHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 

+  NP5,XSHIFT,YSHIFT,NAIRFO,XI(202),YI(202), 

+  COSTHL(201),SINTHL(201) 

COMMON  /PARD/  NACAD(2) ,TAUD(2) ,EPSMAD(2) ,PTMAXD(2) 
READ  (1,501)  ITITLE 
WRITE  (6,501)  ITITLE 
DO  10   I  =  1, ITITLE 
READ  (1,502)  TITLE 
10   WRITE  (6,503)  TITLE 

501  FORMAT(3I5) 

502  FORMAT(20A4) 

503  F0RMAT(1X,20A4) 


100 


504 


F0RMAT(I5,2F10. 6) 
READ  (1,501)  IFLAG,NLOWER,NUPPER 
WRITE  (6,501)  IFLAG,NLOWER,NUPPER 
READ  (1,504)  NAIRFO.XSHIFT, YSHIFT 
WRITE  (6,504)  NAIRFO,XHIFT, YSHIFT 
IF  (I FLAG  .NE.  0)  RETURN 

COMPUTATION  FOR  NACA  4-DIGITS  SERIES 

READ  (1,501)  NACAD(l) 

WRITE  (6,501)  NACAD(l) 

READ  (1,501)  NACAD(2) 

WRITE  (6,501)  NACAD(2) 

DO  100  I=l,NAIRFO 

IEPS      =  NACAD(I)/1000 

IPTMAX    =  NACAD(I)/100  -  10*IEPS 

ITAU      =  NACAD(I)  -  1000*IEPS 

EPSMAD(I)  =  IEPS---0.  01 

PTMAXD(I)   =  IPTMAX*0.  1 

TAUD(I)    =  ITAU* 0.01 

IF  (IEPS  .LT.  10)  GOTO  100 

COMPUTATION  FOR  NACA  5 -DIGITS  SERIES  NOTING  THAT  TAUD(I) 
COMPUTED  AS  PER  4-DIGITS  SERIES. 


100*IPTMAX 


IS 


PTMAXD(I)   =  0. 2025 
EPSMAD(I)  =  2.  6595*PTMAXD(I)**3 
100   CONTINUE 
RETURN 
END 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


SUBROUTINE  SETUP 


SETUP  COORDINATES  OF  PANEL  NODES  AND  SLOPES  OF  PANELS 
COORDINATES  ARE  READ  FROM  INPUT  DATA  FILE  UNLESS 

THE  AIRFOIL  IS  OF  NACA  XXXX  OR  NACA  230XX  TYPE 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C     NACA(I)    DUMMY  NACA  NUMBER  FOR  TRANSFER  TO  SUBROUTINES 

C  BODY  AND  NACA 

C     TAU(I)    DUMMY  MAX  THICKNESS  FOR  TRANSFER  TO  SUBROUTINES 

C  BODY  AND  NACA 

C     EPSMAX(I)  DUMMY  MAX  CAMBER  FOR  TRANSFER  TO  SUBROUTINES 

C  BODY  AND  NACA 

C     PTMAX(I)   DUMMY  CHORDWISE  POSN  FOR  TRANSFER  TO  SUBROUTINES 

C  BODY  AND  NACA 

C     

SUBROUTINE  SETUP 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X( 202) ,Y( 202) , 
+        COSTHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 
+        NP5, XSHIFT, YSHIFT, NAIRFO,XI( 202) , YI( 202) , 
+        COSTHL(201),SINTHL(201) 

COMMON  /PARD/  NACAD( 2) ,TAUD(2) ,EPSMAD(2) ,PTMAXD( 2) 

COMMON  /PAR/  NACA,TAU,EPSMAX,PTMAX 
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COMMON  /NUM/  PI.PI2INV 

PI      =  3. 1415926585 

PI2INV  =  .  5/PI 
C 

C     SET  COORDINATES  OF  NODES  ON  BODY  SURFACE 
C 

DO  210  I=l,NAIRFO 

NODTOT  =  NLOWER  +  NUPPER 

K=(I-l)*(NODTOT+l) 

IF  (IFLAG  .NE.  0)  GO  TO  10 
C 

C     NACA  SERIES  AIRFOIL  CALCULATIONS 
C 

NACA=NACAD(I) 

TAU=TAUD(I) 

EPSMAX=EPSMAD(I) 

PTMAX=PTMAXD(I) 

NPOINT  =  NLOWER 

SIGN    =  -1.0 

NSTART  =  0 

DO  110  NSURF  =1,2 

DO  100  N  =  1, NPOINT 

FRACT   =  FLOAT(N-l)/FLOAT( NPOINT) 

Z       =  . 5*(1.  -  COS(PI*FRACT)) 

J       =  NSTART  +  N 

CALL  BODY(Z,SIGN,XD,YD) 

XI(J+K)  =  XD 

YI(J+K)  =  YD 
100   CONTINUE 

NPOINT  =  NUPPER 

SIGN    =1.0 

NSTART  =  NLOWER 
110  CONTINUE 

XI(I-(NODTOT+l))  =  XI(1+K) 
120  YI(I*(NODTOT+l))  =  YI(1+K) 

GO  TO  210 
C 

C     AIRFOIL  DATA  INPUT  THAT  ARE  NOT  NACA  SERIES  AIRFOIL. 
C 

10   READ  (1,501)  (XI(J+K),J=l,NODTOT+l) 

READ  (1,501)  (YI(J+K),J=l,NODTOT+l) 

WRITE  (6,501)  (XI(J+K),J=l,NODTOT+l) 

WRITE  (6,501)  (YI(J+K),J=l,NODTOT+l) 
501  FORMAT  (6F10. 6) 
210  CONTINUE 

NP1     =  NODTOT  +  1 

NP2  =  NAIRF0*NP1+1 

C 

C  SET  SLOPES  OF  PANELS  AND  CALCULATE  AIRFOIL  PERIMETER 

C 

DO   220      I   =   l,NAIRFO 

K  =  (I-l)*(NODTOT+l) 

SS(I)        =  0.0 

DO  200     J  =   1, NODTOT 

DX  =  XKJ+l+K)    -  XI(J+K) 

DY  =  YKJ+l+K)    -   YI(J+K) 
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DIST    =  SQRT(DX*DX  +DY*DY) 

SS(I)   =  SS(I)  +  DIST 

SINTHL(J+K)  =  DY/DIST 

COSTHL(J+K)  =  DX/DIST 
200  CONTINUE 
220  CONTINUE 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C     SUBROUTINE  BODY( Z ,SIGN,X, Y)  C 

c  c 

C  RETURN  COORDINATES  OF  POINT  ON  THE  BODY  SURFACE  C 

C  C 

C  Z  =  NODE -SPACING  PARAMETER  C 

C  X,Y  =  CARTESIAN  COORDINATES  C 

C  SIGN  =  +1.  FOR  UPPER  SURFACE  C 

C  -1.  FOR  LOWER  SURFACE  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE  BODY( Z ,SIGN,XD,YD) 

COMMON  /PAR/  NACA,TAU,EPSMAX,PTMAX 

IF  (SIGN  .LT.  0.  0)    Z  =  1.  -  Z 

CALL  NACA45(Z, THICK, CAMBER, BETA) 

XD      =  Z  -  SIGN*THICK*SIN(BETA) 

YD      =  CAMBER  +  SIGN*THICK*COS(BETA) 

RETURN 

END 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C     SUBROUTINE  NACA45( Z , THICK, CAMBER, BETA)  C 

C  C 

C  EVALUATE  THICKNESS  AND  CAMBER  C 

C  FOR  NACA  4-  OR  5 -DIGIT  AIRFOIL  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C     Z      ABSCISSA  OF  CAMBERLINE  POINT 

C     X,Y    LOCAL  CARTESIAN  COORDINATES  OF  NACA-AIRFOIL 

C     SIGN    SURFACE  INDICATOR:  SIGN  =  1,  UPPER  SURFACE 

C  =-1,  LOWER  SURFACE 

C     THICK   THICKNESS  AT  THE  POSITION  Z 

C     CAMBER  CAMBER  AT  THE  POSITION  Z 

C     BETA   ANGLE  BETWEEN  CAMBER  LINE  AND  X-AXIS  AT  POSITION  Z 

C     TAU    MAXIMUM  THICKNESS  (INPUT) 

C     EPSMAX  MAXIMUM  CAMBER  ( INPUT) 

C     PTMAX   COORDINATE  POSITION  OF  MAXIMUM  CAMBER  (INPUT) 

C 

SUBROUTINE  NACA45( Z , THICK, CAMBER, BETA) 

COMMON  /PAR/  NACA, TAU, EPSMAX, PTMAX 

THICK   =0.0 

IF  (Z  .LT.  l.E-10)     GO  TO  100 

THICK   =  5.*TAU*(.  2969*SQRT(Z)  -  Z*(.126  +  Z*(.3537 
+  -  Z*(.2843  -  Z*.  1015)))) 

100   IF  (EPSMAX  . EQ.  0.0)  GO  TO  130 
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IF  (NACA  .GT.  9999)  GO  TO  140 
C 

C     CAMBERLINE  OF  NACA  4-DIGIT  SERIES 
C 

IF  (Z  .GT.  PTMAX)  GO  TO  110 
C 

C     FORWARD  PART  OF  CAMBER  LINE 
C 

CAMBER  =  EPSMAX/PTMAX/PTMAX*(2.*PTMAX  -  Z)*Z 

DCAMDX  =  2.*EPSMAX/PTMAX/PTMAX*(PTMAX  -  Z) 

GO  TO  120 
C 

C     AFT  PART  OF  CAMBERLINE 
C 
110  CAMBER  =  EPSMAX/(1.  -PTMAX)**2*( l.  +  Z  -  2.  *PTMAX)*(  1.  -  Z) 

DCAMDX  =  2.*EPSMAX/(1. -PTMAX)** 2* (PTMAX  -  Z) 
120   BETA    =  ATAN( DCAMDX) 

RETURN 
130  CAMBER  =  0.  0 

BETA    =  0.  0 

RETURN 
C 

C     CAMBERLINE  OF  NACA  5 -DIGIT  SERIES 
C 

140   IF  (Z  .GT.  PTMAX)      GO  TO  150 
C 

C     FORWARD  PART  OF  CAMBER  LINE 
C 

W       =  Z/ PTMAX 

CAMBER  =  EPSMAX*W*((W  -  3. )*W  +  3.  -  PTMAX) 

DCAMDX  =  EPSMAX*3.*W*(1.  -  W) /PTMAX 

GO  TO  120 
C 

C     AFT  PART  OF  CAMBERLINE 
C 
150  CAMBER  =  EPSMAX-(1.  -  Z) 

DCAMDX  =  -  EPSMAX 

GO  TO  120 

END 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C     SUBROUTINE  GAUSS(NRHS ,M,NITR)  C 

C  C 

C  SOLUTION  OF  LINEAR  ALGEBRAIC  SYSTEM  BY  C 

C  GAUSS  ELIMINATION  WITHOUT  PARTIAL  PIVOTING  C 

C  C 

C  °A  COEFFICIENT  MATRIX  C 

C  NEQNS     =  NUMBER  OF  EQUATIONS  C 

C  NRHS      =  NUMBER  OF  RIGHT  HAND  SIDES  C 

C  C 

C  RIGHT-HAND  SIDES  AND  SOLUTIONS  STORED  IN  C 

C  COLUMNS  NEQNS+1  THRU  NEQNS+NRHS  OF  °A  C 

c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  GAUSS(NRHS ,M,NITR) 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y(202) , 
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+  C0STHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 

+  NP5,XSHIFT,YSHIFT,NAIRF0,XI(202),YI(202), 

+  COSTHL(201),SINTHL(201) 

COMMON  /COF/  A( 20 1,211), KEQNS 
C     WRITE  (9,900)  M,NITR 

C900  FORMAT  (1X,'M  =',13,'     NITR  =',I3/) 

C     WRITE  (9,910)(I,A(I,NP3),A(I,NP4),A(I,NP5),I=1,2*N0DT0T) 
910   FORMAT  (1X,I5,3E14.  6) 

IF  (M. EQ. 0)   KEQNS  =  NODTOT 
NEQNS   =  KEQNS  *  NAIRFO 
NP      =  NEQNS  +  1 
NTOT    =  NEQNS  +  NRHS 
IF  (NITR  .GT.  0)  GO  TO  160 
C 

C  GAUSS  REDUCTION 

C 

DO  150   I  =  2, NEQNS 
IM      =1-1 
C 

C  ELIMINATE  (I-l)TH  UNKNOWN  FROM 

C  ITH  THRU  ( NEQNS )TH  EQUATIONS 

C 

DO  150  J  =  I, NEQNS 
R       =  A(J,IM)/A(IM,IM) 
DO  150   K  =  I, NTOT 
150  A(J,K)   =  A(J,K)  -  R*A(IM,K) 
GO  TO  170 
C 

C  GAUSSIAN  ELIMINATION  ON  ONLY  THE  RIGHT-HAND-SIDES 

C 
160   DO  180   I  =  2, NEQNS 
IM      =1-1 
DO  180   J  =  I, NEQNS 
R       =  A(J,IM)/A(IM,IM) 
DO  180   K  =  NP.NTOT 
180   A(J,K)   =  A(J,K)  -  R-A(IM,K) 
170  CONTINUE 
C 

C  BACK  SUBSTITUTION 

C 

DO  220   K  =  NP.NTOT 

A( NEQNS, K)  =  A(NEQNS,K)/A(NEQNS, NEQNS) 
DO  210   L  =  2, NEQNS 
I       =  NEQNS  +  1  -  L 
IP      =1+1 
DO  200  J  =  IP, NEQNS 
200  A(I,K)   =  A(I,K)  -  A(I,J)*A(J,K) 
210  A(I,K)   =  A(I,K)/A(I,I) 
220  CONTINUE 
RETURN 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C     SUBROUTINE  VELDIS(SINALF,COSALF)  C 

C  C 

C     COMPUTE  STEADY  FLOW  PRESSURE  DISTRIBUTION  AND  VELOCITY  C 
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C     POTENTIAL  AT  MID-POINTS  OF  PANELS  FOR  THE  STEADY  FLOW  CASE       C 
C  C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C     VTANG  TANGENTIAL  VELO  COMP  OF  A  FIXED  POINT  OF  THE  MOVING 

C  LOCAL  FRAME  OF  REFERENCE. 

C     PHI(I)  DISTURBANCE  VELO  POTENTIAL  AT  THE  MID  POINT  OF  THE 

C  I-TH  PANEL 

C     PHILE(L) DIFFERENCE  OF  THE  POTENTIALS  OF  THE  LEADING  EDGE  TO 

C  THE  LOWER  TRAILING  EDGE  FOR  THE  RESPECTIVE  AIRFOIL 

C     PIN(L)   DIFFERENCE  OF  THE  POTENTIALS  AT  A  POINT  1000  CHORD 

C  LENGTH  UPSTREAM  OF  THE  LE  FOR  THE  RESPECTIVE  AIRFOIL 

C     SUMC(L)  GAMMA  ASSOCIATED  WITH  THE  INTEGRATION  OF  THE  DISTURB- 

C  ANCE  VELOCITY  AROUND  THE  WHOLE  AIRFOIL 

C     

SUBROUTINE  VELDIS 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y(202) , 
+  COSTHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 

+  NP5,XSHIFT,YSHIFT,NAIRFO,XI(202),YI(202), 

+  COSTHL(201),SINTHL(201) 

COMMON  /COF/  A(201,211),KEQNS 
COMMON  /CPD/  CP(200),SCL,T,SCM,SGAM 
COMMON  /NUM/  PI,PI2INV 

COMMON  /SING/  Q( 200) ,GAMMA( 2) ,QK(200) ,GAMK(2) 
COMMON  /POT/  PHI(200),PHIK(200) 
COMMON  /GUST/  UG( 200) , VG( 200) ,XGF ,UGUST, VGUST 
COMMON  /EXTV/  UE( 200) , VN( 200) 

COMMON  /GEOM/  SINALF( 2) ,COSALF( 2) ,OMEGA( 2) ,UX( 2) ,UY( 2) , PIVOT, 
+  XPRM,YPRM 

DIMENSION  CKUTTA(2),PHITEL(2),PHITEU(2) 
7770  DIMENSION  CONTR2(2) 
6666  REAL  *8  PIN(2),VELX 

DIMENSION  WGHT(5) ,PLOC(5) ,SUMC( 2) , 
+AANPK 50, 50,5 ),AANP2( 50,50,5 ) ,BBNP1( 50,50,5) ,BBNP2( 50, 50,5), 
+COSTHP( 102,6) ,SINTHP( 102,6) ,AANP4(2) ,PHIL( 102) ,UGU( 100,6) , 
+VGU(100,6),PHILE(2) 

LOCATION  AND  WEIGHTING  VALVES  FOR  THE  GAUSSIAN  QUADRATURE  USED 
TO  INTEGRATE  FOR  THE  VELO  POTENTIAL  AROUND  THE  AIRFOIL 

DATA  WGHT/. 11846344, . 23931434, . 28444444, 
+  . 23931434,. 11846344/ 

DATA  PLOC/. 04691008,.  23076535,.  50000000, 
+  .76923466,-95300899/ 

TH  PANEL 


c 
c 
c 

FIND  VT 

AND  CP  AT  MID-POINT  OF  I- 

DO  140 

L=   l,NAIRFO 

KI 

=  (L-1)*NP1 

LI 

=  (L-l)*NODTOT 

SUMC(L) 

=  0.0 

DO  130 

I  =  l,NODTOT 

XMID 

=  . 5*(X(I+KI)  +X(I+KI+1)) 

YMID 

=  .5*(Y(I+KI)  +Y(I+KI+1)) 

DX 

=  X(I+KI+1)-X(I+KI) 

DY 

=  Y(I+KI+1)-Y(I+KI) 
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100 


120 
150 


155 


130 
140 


DIST  =  SQRT(DX*DX+DY*DY) 

VTANG  =  CCSALFCLV'COSTHLCI+KI) 

:R>:  =  SINALF(L)*COSTHL(I+KI) 

VTFREE  =  VTANG 

VACT  =  VTANG 


SINALF(L)*SINTHL(I+KI) 
COSALF(L)*SINTHL(I+KI) 


ADD  CONTRIBUTION  OF  J-TH  PANEL 

1.  CONTRIBUTION  OF  J-TH  PANEL  TO  THE  VELO  COMP  OF  THE  MIDPT 
OF  THE  I-TH  PANEL 

2.  CONTRIBUTION  TO  THE  VELO  POTENTIAL.  THIS  IS  DONE  BY 
INTEGRATING  OVER  SMALLER  PANELS  OF  THE  AIRFOIL. 


::: 


DO  155  K 

DX 

DY 

XMID 

ym:: 

VDUM 

DO  150 

KJ 

LJ 

DO  120 

FLOG 

FTAN 

IF  (J+KJ 

DXJ 

■ 
DYJ 
DYJP 
FLOG 
FTAN 

ctimtj  : 

st:::t: 

AA 

B 

VDUM 

IF(K 

IF(K 

CONTINUE 

CONTINUE 

VTANG 

SUMC(L)  = 

CONTINUE 

PHKI+LI) 

CP(I+LI)= 

UE(I+LI)= 

VN(I+LI)= 

CONTINUE 

CONTINUE 


■X(I+KT)) 
■Y(I+KI)) 


EQ. 

EQ. 


=  1,5 

PLOC(K)*(X(I+KI+l) 

PLOC(K)*(Y(I+KI+l) 

X(I+KI)  +  DX 
Y(I+KI)  +  DY 
0.0 
=  l.NAIRFO 
CLM-1)*NP1 
(LM-1)*NODTOT 
=  l.NODTOT 
0.  0 

PI 

EQ.  I+KI)  GO  TO  100 

XMID  -  X(J+KJ) 

XMID  -  X(J+KJ+1) 

YMID  -  Y(J+KJ) 

YMID  -  Y(J+KJ+1) 

.  5*ALOG(  ( DXJP*DXJP+DYJP*DYJP) / ( DXJ*DXJ+DYJ*DYJ) ) 

ATAN2(  DYJP-'VDXJ-DXJP^DYJ ,  DXJP*DXJ+DYJP*DYJ) 

COSTHEC  I+KI  )--'--C0STHE(  J+KJ)  +  SINTHE(  I+KI  )*SINTHE(  J+KJ) 

SINTHEC I+KI )*COSTHE( J+KJ)  -  COSTHEC I+KI )*SINTHE( J+KJ) 

PI2INV*(FTAN*CTIMTJ  +  FLOG'vSTIMTJ) 

PI2INV*(FXOG*CTIMTJ  -  FTAN* ST I MT J) 

VDUM  -  B*Q(J+LJ)  +GAMMA(LM)*AA 

3)  VACT  =  VACT  -  B*Q(J+LJ)  +  GAMMA(LM)*AA 

3)  VNORM  =  VNORM  +  AA*Q(J+LJ)  +  GAMMA (LM)*B 


VTANG  +VDUM  *VGHT(K) 

:(L)  +VDUM*DIST*WGHT(K) 


=  (VTANG-\TFREE)*DIST 
1.0  -  VACT*VACT 
VACT 
VNORM 


COMPUTE  DISTURBANCE  POTENTIAL  AT  THE  LEADING  EDGE  BY  LINE 
INTEGRAL  OF  THE  VELOCITY  FIELD 
FROM  UPSTREAM  (AT  INFINITY)  TO  THE  LEADING  EDGE 


:  n 


L  =  l.NAIRFO 
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20 

40 

30 
55 


YMID  =  PIVOT  *  SINALF(L)  +  (L-1)*YSHIFT 

XLE  =  -PIV0T*C0SALF(L)+(L-1)*XSHIFT 

XL  =  XLE 

XL  =  0.  0 

NPHI  =  10  *  NLOWER 

PIN(L)   =0.0 

DO  30  I  =  l.NPHI 

FRACT  =  FLOAT(I)/FLOAT(NPHI) 

XLP  =  -100.0  *  (1.0  -  COS(0.  5*PI*FRACT))+XLE 

IF  (I  .EQ.  1)  XLP  =  -0.000197+XLE 

XLP  =  -10.0  *  (1.0  -  C0S(0.5*PI*FRACT)) 

DELX  =  XL  -  XLP 

XMID  =  0. 5*(XL+XLP) 

XMID  =  0.5*(XL+XLP)*COSALF(L) 

YMID  =  0. 5*(XL+XLP)*SINALF(L) 

XL  =  XLP 

VELX  =  UGUST 

ADD  CONTRIBUTION  OF  J-TH  PANEL 


DO  40 
KJ 

LJ 

DO  20 

DXJ 

DXJP 

DYJ 

DYJP 

FLOG 

FTAN 

CALMTJ 

SALMTJ 

APY 

BPY 

VELX 

CONTINUE 

CONTINUE 

PIN(L) 

CONTINUE 

CONTINUE 


LN  =  l,NAIRFO 

=  (LN-1)*NP1 

=  (LN-l)*NODTOT 

J  =  l,NODTOT 

=  XMID  -  X(J+KJ) 

=  XMID  -  X(J+KJ+1) 

=  YMID  -  Y(J+KJ) 

=  YMID  -  Y(J+KJ+1) 

=  .  5*ALOG((DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 

=  ATAN2 ( DY JP- DX J -DX JP*DY J , DX JP*DX J+DY JP*DY J ) 

=  -COSALF(L)*COSTHE(J+KJ)  -  SINALF(L)*SINTHE( J+KJ) 

=  -SINALF(L)*COSTHE(J+KJ)  +  COSALF(L)*SINTHE( J+KJ) 

=  PI2INV*(FTAN*CALMTJ  +  FLOG*SALMTJ) 

=  PI2INV*(FLOG*CALMTJ  -  FTAN*SALMTJ) 

=  VELX  -  DPROD(BPY,Q(J+LJ))  +DPROD ( GAMMA ( LN ), APY) 


=  PIN(L)  +  VELX  *  DBLE(DELX) 


230 


COMPUTATION  OF  THE  VELOCITY  POTENTIAL  FOR  MIDPOINT  OF  EACH  PANEL 

DO  240  L   =  l,NAIRFO 

LI        =  (L-l)*NODTOT 

PHP        =  -PIN(L) 

BEGIN  WITH  LOWER  SURFACE 

DO  230  I   =  NLOWER, 1,-1 

PHC  =  PHP-PHI (I+LI) 

PHKI+LI)  =  0.5*(PHP+PHC) 

PHP  =  PHC 

PHITEL(L)  =  PHC 
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C     RESET  FOR  UPPER  SURFACE 
C 

PHP       =  -PIN(L) 

DO  250  I   =  NLOWER+1, NODTOT 

PHC       =  PHP+PHKI+LI) 

PHKI+LI)  =  0.5*(PHP+PHC) 
250  PHP       =  PHC 

PHITEU(L)  =  PHC 
240  CONTINUE 

DO  334  L  =  1  ,  NAIRFO 

WRITE  (6,1010)  L 

WRITE  (6,1000) 

KI       =  (L-l)  *  NP1 

LI       =  (L-l)  *  NODTOT 

SUMC(L)   =  SUMC(L)/SS(L) 

DO  333  I  =  1  ,  NODTOT 

XMID    =  .5*(X(I+KI)  +  X(I+KI+1)) 

YMID    =  .5*(Y(I+KI)  +  Y(I+KI+1)) 

333  WRITE  (6,1050)  I+LI , XMID, YMID, Q( I+LI) ,GAMMA(L) ,CP( I+LI) ,UE( I+LI) , 
+  VN( I+LI ),PHI( I+LI ),SUMC(L) 

334  WRITE  (6,235)  PIN(L) 

235  FORMAT  (IX,1 PHI  AT  LEADING  EDGE  =  '  F10.6,/) 

1000  F0RMAT(/,4X,'J* ,4X, 'X( J) ' ,5X. ' Y( J) f,6X, '  Q(  J) ' ,6X, 'GAMMA' ,5X, 

+  'CP(J)' ,6X,'V(J)',6X,'VN(J)f ,6X,*PHI',4X,'INTGAMMA' ,/) 
1010  FORMAT(//'  AIRFOIL  NUMBER' , 14,/) 
1050  FORMAT(I5,9F10. 6) 
RETURN 
END 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C     SUBROUTINE  FANDM(SINALF,COSALF)  C 

C  C 

C  COMPUTE  AND  PRINT  OUT  CD,CL,CM  C 

C  INTEGRATE  PRESSURE  DISTRIBUTION  BY  TRAPEZOIDAL  RULE      C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C     CP(I)  PRESSURE  COEFFICIENT  OF  THE  I-TH  PANEL 

C     CL(L)  COEFFICIENT  OF  LIFT  FOR  THE  L-TH  AIRFOIL 

C     CD(L)  COEFFICIENT  OF  DRAG  FOR  THE  L-TH  AIRFOIL 

C     CM(L)  COEFFICIENT  OF  MOMENT  FOR  THE  L-TH  AIRFOIL  WITH 

C  RESPECT  TO  THE  LEADING  EDGE 

C     CFX(L)  COEFFICEIENT  OF  TOTAL  FORCE  IN  X-DIR  OF  GLOBAL  SYS. 

C     CFY(L)  COEFFICEIENT  OF  TOTAL  FORCE  IN  Y-DIR  OF  GLOBAL  SYS. 

C     

C 

SUBROUTINE  FANDM( SIN1 ,SIN2 ,C0S1 ,COS2) 

COMMON  /BOD/  IFLAG, NLOWER,NUPPER, NODTOT, X( 202 ) ,Y( 202) , 
+  C0STHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 

+  NP5, XSHIFT,YSHIFT, NAIRFO, XI(202) ,YI(202) , 

+  COSTHL(201),SINTHL(201) 

COMMON  /CORV/  CV(2,200)  ,XC(2,200)  ,YC(2,200)  ,M,TD,CWX(2,200) , 
+  CWY(2,200),XCI(2,200),YCI(2,200) 

COMMON  /CPD/  CP(200),SCL,T,SCM,SGAM 
COMMON  /SING/  Q( 200) ,GAMMA( 2) ,QK(200) ,GAMK(2) 
COMMON  /GEOM/  SINALF(2) ,COSALF(2) ,OMEGA( 2) ,UX(2) ,UY( 2) , PIVOT, 
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+  XPRM , YPRM 

DIMENSION  CM(2),CD(2),CL(2),CFX(2),CFY(2) 
C 

C     INITIALISE  COEFFICIENT 
C 

DO  110  L  =  l,NAIRFO 
CM(L)   =0.0 
CFX(L)   =0.0 
110  CFY(L)   =0.0 
C 

C     INTEGRATE  AROUND  THE  AIRFOIL  TO  GET  GLOBAL  X  AND  Y  FORCES 
C 

DO  120  L  =  l.NAIRFO 
KI      =  (L-l)*(NODTOT+l) 
LI  =  (L-l)*NODTOT 
DO  100   I  =  l,NODTOT 
XMID    =  . 5*(XI(I+KI)  +  XId+KI+1)) 
YMID    =  .5*(YI(I+KI)  +  YIU+KI+1)) 
DX      =  XId+KI+1)  -  XI(I+KI) 
DY      =  YKI+KI+1)  -  YI(I+KI) 
CFX(L)   =  CFX(L)  +  CP(I+LI)*DY 
CFY(L)   =  CFY(L)  -  CP(I+LI)*DX 
CM(L)   =  CM(L)  +  CP(I+LI)*(DX*XMID  +  DY*YMID) 
100  CONTINUE 
120   CONTINUE 
C 

C     DECOMPOSE  INTO  LIFT  AND  DRAG  COMPONENTS  W. R. T.  RESP  LOCAL  SYSTEM 
C 

DO  130  L  =l,NAIRFO 

CD(L)   =  CFX(L)*COSALF(L)  +  CFY(L)*SINALF(L) 
CL(L)   =  CFY(L)*COSALF(L)  -  CFX( L)*SINALF(L) 
130  WRITE  (6,1000)  L,CD(L) ,CL(L) ,CM(L) 
IF  (M  .EQ.  0)  RETURN 
IF  (SCL  .NE.  0.  0)  THEN 

WRITE  (8,1100)  T,CL(1)/SCL,CL(2)/SCL,CM(1)/SCM,CM(2)/SCM,CD(1), 
+CD(2),GAMK(1)/SGAM,GAMK(2)/SGAM 
ELSE 

WRITE  (8,1100)  T,CL(1),CL(2),CM(1),CM(2),CD(1),CD(2) 
END  IF 
1000  FORMAT(//,'  AIRFOIL  NO  ',14,/,'     CD  =',F10.6, 

+  '     CL=',F10.  6,'     CM=',F10.6) 

1100  FORMAT(9F10. 6) 
RETURN 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C     SUBROUTINE  INFL  (NITR)  C 

C  C 

C  CALCULATE  INFLUENCE  COEFFICIENTS  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  INFLUENCE  COEFFICIENTS  ON  THE  AIRFOIL  DUE  TO  THE  AIRFOIL: 

C  AAN(I,J)   NORMAL  VELO  AT  THE  MIDPOINT  OF  THE  I-TH  PANEL  DUE 

C  TO  A  SOURCE -D I ST  OF  UNIT  STRENGTH  ON  THE  J-TH  PANEL 

C  SUMAAN(I,L)  ..  NORMAL  VELO  AT  THE  MIDPOINT  OF  THE  I-TH  PANEL  DUE 
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c 

c 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


BBN(I,J)   ... 
SUMAAN(I,L)  . 


SOURCE  DIST  OF  UNIT  STRENGTH  FROM  THE  L  AIRFOIL 
NORMAL  VELO  AT  THE  MIDPOINT  OF  THE  I-TH  PANEL  DUE 
TO  A  VORTEX -DIST  OF  UNIT  STRENGTH  ON  THE  J-TH  PANEL 
NORMAL  VELO  AT  THE  MIDPOINT  OF  THE  I-TH  PANEL  DUE 
VORTEX  DIST  OF  UNIT  STRENGTH  FROM  THE  L  AIRFOIL 


INFLUENCE  COEFFICIENT  ON  THE  WAKE  ELEMENT: 


AYNP1(L,J) 


AYNP1(L,NP3) 


BYNP1(L,J) 


BYNP1(L,NP3) 


CYNP1(L,N) 


CXNP1(L,N) 


Y-VELO  COMP  AT  THE  MIDPOINT  OF  THE  WAKE  PANEL  FROM 
THE  L-TH  AIRFOIL  DUE  TO  A  SOURCE  DIST  OF  UNIT  STRE- 
NGTH FROM  THE  J-TH  PANEL 

Y-VELO  COMP  AT  THE  MIDPOINT  OF  THE  WAKE  PANEL  FROM 
THE  L-TH  PANEL  DUE  TO  A  SOURCE  DIST  OF  UNIT  STREN- 
GTH FROM  THE  WAKE  PANEL  OF  THE  OTHER  AIRFOIL. 
(USED  ONLY  FOR  BXNP1(L,NP3)  SINCE  THERE  IS  NO  SOURC 
DIST  ON  THE  WAKE  PANEL) 

Y-VELO  COMP  AT  THE  MIDPOINT  OF  THE  WAKE  PANEL  FROM 
THE  L-TH  AIRFOIL  DUE  TO  A  VORTEX  DIST  OF  UNIT  STRE- 
NGTH FROM  THE  J-TH  PANEL 

Y-VELO  COMP  AT  THE  MIDPOINT  OF  THE  WAKE  PANEL  FROM 
THE  L-TH  AIRFOIL  DUE  TO  A  VORTEX  DIST  OF  UNIT  STR- 
ENGTH FROM  THE  WAKE  PANEL  OF  THE  OTHER  AIRFOIL. 
Y-VELO  COMP  AT  THE  MIDPOINT  OF  THE  WAKE  PANEL  FROM 
THE  L-TH  AIRFOIL  DUE  TO  THE  N-TH  CORE  VORTEX  OF 
UNIT  STRENGTH 

X-VELO  COMP  AT  THE  MIDPOINT  OF  THE  WAKE  PANEL  FROM 
THE  L-TH  AIRFOIL  DUE  TO  THE  N-TH  CORE  VORTEX  OF 
UNIT  STRENGTH 


C 

c 
c 
c 


INFLUENCE  COEFFICIENTS  ON  THE  AIRFOIL  DUE  TO  THE  WAKE 
SUMCCN(I)  ...  NORMAL  VELO  AT  THE  MIDPOINT  OF  THE  I-TH  PANEL  DUE 

TO  ALL  POINT  VORTICES  OF  ACTUAL  STRENGTH. 
SUMCCT(I)  ...  TANGENTIAL  VELO  AT  THE  MIDPOINT  OF  THE  I-TH  PANEL 

DUE  TO  ALL  POINT  VORTICES  OF  ACTUAL  STRENGTH. 


SUBROUTINE  INFL  (NITR) 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y(202) , 
+  COSTHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 

+  NP5,XSHIFT,YSHIFT,NAIRFO,XI(202)JYI(202), 

+  COSTHL(201),SINTHL(201) 

COMMON  /NUM/  PI,PI2INV 

COMMON  /WAK/  VYW( 2) ,VXW( 2) ,WAKE( 2) ,DT 

COMMON  /CORV/  CV( 2,200) ,XC( 2 ,200) ,YC(2 ,200) ,M,TD,CVVX( 2,200) , 
+  CWY(2,200),XCI(2,200),YCI(2,200) 

COMMON  /INF1/  AAN(201,201),BBN(201,201),AYNP1(2,201),BYNP1(2,201), 
+  SUMAAN(201,2),SUMBBN(201,2) 

COMMON  /INF2/  SUMCCN(201) ,SUMCCT(201) ,CYNP1(2,200) , 
+  CXNP1(2,200) 

COMMON  /GEOM/  SINALF(2) ,COSALF(2) ,OMEGA(2) ,UX(2) ,UY(2) , PIVOT, 
+  XPRM,YPRM 

DIMENSION  GAYNP1(2,201),GBYNP1(2,201) 

IF  (M  .GT.  1)   GO  TO  510 

INFLUENCE  COEFFICIENT  ON  THE  I-TH  PANEL  BY  THE  J-TH  PANEL  FROM 

THE  SAME  AIRFOIL. 


Ill 


100 


110 
120 
200 
510 


DO  200  L  =  l,NAIRFO 

LI      =  (L-l)*NODTOT 

KI      =  (L-1)*NP1 

DO   120      I   =   l,NODTOT 

XMID  =  .  5*(XI(I+KI)   +  XKI+KI+l)) 

YMID  =  .5*(YI(I+KI)   +  YKI+KI+1)) 

SUMAAN(I+LI,L)   =  0.  0 

SUMBBN(I+LI,L)   =0.0 

DO  110  J  =  l,NODTOT 

0.  0 

PI 

EQ.  I+LI)     GO  TO  100 

XMID  -  XKJ+KI) 

XMID  -  XKJ+KI+1) 

YMID  -  YICJ+KI) 

YMID  -  YI(J+KI+1) 

. 5  * ALOG ( ( DX JP*DX JP+D Y JP*D Y JP ) / ( DX J*DX J+D Y J*D Y J ) ) 

ATAN2 ( D Y JP*DX J - DX JP*D Y J , DX JP*DX J+D Y JP*D Y J ) 

COSTHL(I+KI)*COSTHL(J+KI)  +  SINTHL(  I+KI)*SINTHL(  J+KI ) 

SINTHL(I+KI)*COSTHL(J+KI)  -  COSTHL( I+KI)*SINTHL( J+KI) 
PI2INV*(FTAN*CTIMTJ  +  FLOG*STIMTJ) 
PI2INV*(FLOG*CTIMTJ  -  FTAN*STIMTJ) 
SUMAAN(I+LI,L)  +  AAN( I+LI , J+LI) 
SUMBBN(I+LI,L)  +  BBN( I+LI , J+LI) 


FLOG 

FT  AN 

IF  (J+LI 

DXJ 

DXJP 

DYJ 

DYJP 

FLOG 

FT  AN 

CTIMTJ 

STIMTJ 

AAN( I+LI, J+LI) 

BBN( I+LI, J+LI) 

SUMAAN(I+LI,L) 

SUMBBN(I+LI,L) 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

IF  (NITR 


GT.  0)  GO  TO  271 


INFLUENCE 
THE  OTHER 


270 


COEFFICIENT 
AIRFOIL. 


ON  THE  I-TH  PANEL  BY  THE  J-TH  PANEL  FROM 


DO 

LI 

KI 

DO  260 

XMID 

YMID 

SUMAAN( I+LI, 3-L)  = 

SUMBBN(I+LI,3-L)  = 


=  1,NAIRF0 
(L-1)*N0DT0T 
(L-l)-NPl 
=  1,N0DT0T 
.  5*(X(I+KI)  + 
.  5*(Y(I+KI)  + 
0.0 
0.  0 


XCI+KI+1)) 
Y(I+KI+1)) 


DO  250 

DXJ 

DXJP 

DYJ 

DYJP 

FLOG 

FTAN 

CTIMTJ 

STIMTJ 


=  NODTOT+2, 
XMID  -  X(J- 
XMID  -  X(J- 
YMID  -  Y(J- 
YMID  -  Y(J- 


2*NODTOT+l 

KI) 

KI+1) 

KI) 

KI+1) 


250 
260 


AAN(I+LI,J-LI-1) 
BBN(I+LI,J-LI-1) 
SUMAAN( I+LI, 3-L) 
SUMBBN( I+LI, 3-L) 

CONTINUE 
CONTINUE 


. 5* ALOG ( (DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 
ATAN2( DYJP*DXJ-DXJP*DYJ , DXJP*DXJ+DYJP*DYJ) 
COSTHE(I+KI)*COSTHE(J-KI)  +  SINTHE( I+KI)*SINTHE( J-KI) 
SINTHE(I+KI)*COSTHE(J-KI)  -  COSTHE( I+KI)*SINTHE( J-KI) 
PI2INV*(FTAN*CTIMTJ  +  FLOG*STIMTJ) 
PI2INV*(FLOG*CTIMTJ  -  FTAN*STIMTJ) 
SUMAANC I+LI, 3-L)  +  AAN( I+LI , J-LI-1) 
SUMBBN( I+LI, 3-L)  +  BBN( I+LI, J-LI-1) 
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270  CONTINUE 


271 


130 


END  OF  STEADY  FLOW  CALCULATION  FOR  INFLUENCE  COEFFICIENT. 

IF  (M.EQ.  0)   RETURN 
CONTINUE 

INFLUENCE  COEFFICIENT  ON  THE  WAKE  ELEMENT  FROM  THE  AIRFOIL. 

DO  130  L  =  l,NAIRFO 

I       =  NP1*L 

XMID    =  . 5*(X(I)  +  X(NAIRF0*NP1+L)) 

YMID    =  .5*(Y(I)  +  Y(NAIRF0*NP1+L)) 

DO  130  MM  =  l,NAIRFO 

LJ      =(MM-l)*NODTOT 

KJ      =(MM-1)*NP1 

DO  130  J  =  l.NODTOT 

DXJ     =  XMID  -  X(J+KJ) 

DXJP    =  XMID  -  X(J+KJ+1) 

DYJ     =  YMID  -  Y(J+KJ) 

DYJP    =  YMID  -  Y(J+KJ+1) 

FLOG    =  .  5*ALOG(  ( DXJP*DXJP+DYJP*DYJP)  /( DXJ*DXJ+DYJ*DYJ) ) 

FTAN    =  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 

CTIMTJ  =  COSTHE(I)*COSTHE(J+KJ)  +  SINTHE(I)*SINTHE( J+KJ) 

STIMTJ  =  SINTHE(I)*COSTHE(J+KJ)  -  COSTHE( I)*SINTHE( J+KJ) 

GAYNP1(L,J+LJ)  =  PI2INV*(FTAN*COSTHE(J+KJ)  -  FLOG*SINTHE( J+KJ) ) 

GBYNP1(L,J+LJ)  =  PI2INV*(FLOG*COSTHE(J+KJ)  +  FTAN*SINTHE( J+KJ) ) 

AYNP1(L,J+LJ)  =  GAYNP1(L,J+LJ)*C0SALF(L)-GBYNP1(L,J+LJ)*SINALF(L) 

BYNP1(L,J+LJ)  =  GAYNP1(L,J+LJ)*SINALF(L)+GBYNP1(L,J+LJ)*C0SALF(L) 

CONTINUE 

INFLUENCE  COEFICIENT  OF  WAKE  ELEMENT  DUE  TO  WAKE  ELEMENT  FROM 
THE  OTHER  AIRFOIL 


DO  300   L 

I 

XMID 

YMID 

KJ 

MJ 

NJ 

DXJ 

DXJP 

DYJ 

DYJP 

FLOG 

FTAN 

CTIMTJ  = 

STIMTJ  = 


X(NAIRF0*NP1+L)) 
Y(NAIRF0-NP1+L)) 


300 


GAYNP1(L,NP3) 

GBYNP1(L,NP3) 

AYNP1(L,NP3) 

BYNP1(L,NP3) 

CONTINUE 


=  l.NAIRFO 

NP1*L 

.5*(X(I)  + 

,5*(Y(I)  + 

CL-1)*NP1 

2*NP1 

MJ+3 

XMID  -  X(MJ-KJ) 

XMID  -  X(NJ-L) 

YMID  -  Y(MJ-KJ) 

YMID  -  Y(NJ-L) 

.  5*ALOG(  ( DXJP*DXJP+DYJP*DYJP) / ( DXJ*DXJ+DYJ*DYJ) ) 

AT AN 2 ( DY JP* DX J - DX JP*DY J , DX JP*DX J+DY JP*DY J ) 

COSTHE(I)*COSTHE(MJ-KJ)  +  SINTHE(I)*SINTHE(MJ-KJ) 

SINTHE(I)*COSTHE(MJ-KJ)  -  COSTHE( I)*SINTHE(MJ-KJ) 

PI2INV*(FTAN*COSTHE(MJ-KJ)  -  FLOG»SINTHE(MJ-KJ)) 
PI2INV*(FL0G*C0STHE(MJ-KJ)  +  FTAN»SINTHE(MJ-KJ)) 
GAYNP1(L,NP3)*C0SALF(L)-GBYNP1(L,NP3)*SINALF(L) 
GAYNP1(L,NP3)*SINALF(L)+GBYNP1(L,NP3)*C0SALF(L) 
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INFLUENCE  COEFFICIENT  ON  THE  AIRFOIL  BY  THE  WAKE  ELEMENT. 


145 
140 
160 


DO  160 

LI=(L-1) 

KI=(L-1) 

DO  140 

XMID 

YMID 

DO  145 

J 

DXJ 

DXJP 

DYJ 

DYJP 

FLOG 

FTAN 

CTIMTJ 

STIMTJ 

AAN( I+LI 

BBN(I+LI 

CONTINUE 

CONTINUE 

CONTINUE 


L  =  l,NAIRFO 

*NODTOT 

*NP1 

I  =  l,NODTOT 

=  . 5*(X(I+KI)  + 

=  . 5*(Y(I+KI)  + 

MM  =  1,  NAIRFO 

=  NP1*MM 

=  XMID  - 

=  XMID  - 

=  YMID  - 

=  YMID  - 


X(I+KI+1)) 
Y(I+KI+1)) 


X(J) 

X(NAIRF0*NP1+MM) 

Y(J) 

Y(NAIRF0*NP1+MM) 
=  . 5*ALOG( (DXJP*DXJP+DYJP*DYJP) /(DXJ*DXJ+DYJ*DYJ) ) 
=  ATAN2  ( DY  JP*DX  J  -DX  JP*DY  J ,  DX  JP*  DX  J+DY  JP*DY  J ) 
=  COSTHE(I+KI)*COSTHE(J)  +  SINTHE( I+KI)*SINTHE( J) 
=  SINTHE(I+KI)*COSTHE(J)  -  COSTHE( I+KI)*SINTHE( J) 
,2*NODTOT+MM)  =  PI2INV*(FTAN*CTIMTJ  +  FLOG*STIMTJ) 
,2*NODTOT+MM)  =  PI2INV»(FLOG*CTIMTJ  -  FTAN*STIMTJ) 


IF  (M.  EQ.  1)  RETURN 
MM1     =  M  -  1 

INFLUENCE  COEFFICIENT  ON  THE  WAKE  ELEMENT  BY  THE  CORE  VORTICES. 


DO  350   L  =  1, NAIRFO 


XMID 
YMID 
DO  240 


X(NP1*NAIRF0+L)) 
Y(NP1*NAIRF0+L)) 


=  0.5*(X(NP1*L) 
=  0.5*(Y(NP1*L) 
MM  =  1,  NAIRFO 

KN      =  (MM-1)*MM1 

DO  230  N  =  1,MM1 

DX      =  XMID  -  XC(MM,N) 

DY      =  YMID  -  YC(MM,N) 

DIST2   =  DX*DX+DY*DY 

GCYNP1  =  -PI2INV*DX/DIST2 

GCXNP1  =  +PI2INV--DY/DIST2 

CYNP1(L,N+KN)  =  GCYNP1*C0SALF(L)+GCXNP1*SINALF(L) 

CXNP1(L,N+KN)  =  -GCYNP1*SINALF(L)+GCXNP1*C0SALF(L) 
230  CONTINUE 
240  CONTINUE 
350  CONTINUE 

IF  (NITR.  GT.  0)  RETURN 

3     INFLUENCE  COEFFICIENT  ON  THE  AIRFOIL  BY  THE  CORE  VORTICES 


DO  400  L  =  1, NAIRFO 

LI  =(L-l)*NODTOT 

KI  =(L-1)*NP1 

DO  220  I  =  l,NODTOT 

XMID  =  0.5*(X(I+KI) 

YMID  =  0.5*(Y(I+KI) 


X(I+KI+1)) 
Y(I+KI+1)) 
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SUMCCN(I+LI)  =  0.  0 
SUMCCT(I+LI)  =  0.  0 
DO  280  MM  =  1.  NAIRFO 
KN      =  (MM-1)*MM1 
DO  210  N  =  1,MM1 
DX      =  XMID  -  XC(MM,N) 
DY      =  YMID  -  YC(MM,N) 
DIST    =  SQRT(DX*DX+DY*DY) 
COSTHN  =  DX/DIST 
SINTHN  =  DY/DIST 

CTIMTN  =  COSTHE(I+KI)*COSTHN  +  SINTHE(I+KI)*SINTHN 
STIMTN  =  SINTHE(I+KI)*COSTHN  -  COSTHE(I+KI)*SINTHN 
CCN     =  -CTIMTN/DIST 
CCT     =  -STIMTN/DIST 

SUMCCNCI+LI)  =  SUMCCN(I+LI)  +  CCN*CV(MM,N) 
SUMCCT(I+LI)  =  SUMCCT(I+LI)  +  CCT*CV(MM,N) 
CONTINUE 
CONTINUE 

SUMCCNCI+LI)  =  PI2INV*SUMCCN(I+LI) 
SUMCCT(I+LI)  =  PI2INV*SUMCCT(I+LI) 
CONTINUE 
CONTINUE 
RETURN 
END 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


210 
280 


220 
400 


SUBROUTINE  COEF  (SINALF,COSALF, OMEGA, UX,UY,NITR) 

SET  COEFFICIENTS  OF  N  EQUS  ARISING  FROM  FLOW 

TANGENCY  CONDITIONS  AT  MID  POINTS  OF  PANELS 

SOLVING  THE  N- SOURCE  STRENGTHS  IN  TERMS  OF  THE 
VORTICITY  STRENGTH  (RESULTING  IN  2  RHS) 

KUTTA  CONDITION  IS  SATISFIED  SEPARATELY  TO  OBTAIN 
THE  VORTICITY  STRENGTH 

THIS  SOLUTION  METHOD  IS  DESIRED  FOR  UNSTEADY  FLOW 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C  A(I,J)   LHS,  NORMAL  VELOCITY  AT  THE  MIDPOINT  OF  THE  I-TH 

C  PANEL  INDUCED  BY  UNIT  SOURCE  DIST  ON  THE  J-TH  PANEL 

C  A(I,NP3)  FIRST  RHS,  COMPONENT  OF  NORMAL  VELO  AT  THE  I-TH 

C  PANEL  WHICH  IS  DEPENDENT  ON  VORTICITY  STRENGTH  OF 

C  THE  FIRST  AIRFOIL 

C  A(I,NP4)  SECOND  RHS,  COMPONENT  OF  NORMAL  VELO  AT  THE  I-TH 

C  PANEL  WHICH  IS  DEPENDENT  ON  VORTICITY  STRENGTH  OF 

C  THE  SECOND  AIRFOIL 

C  A(I,NP5)  THIRD  RHS,  COMPONENT  OF  NORMAL  VELO  AT  THE  I-TH 

C  PANEL  WHICH  IS  INDEPENDENT  OF  THE  CIRCULATION  OF 

C  BOTH  AIRFOILS 

C  

SUBROUTINE  COEF  (NITR) 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y( 202) , 
+  COSTHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 

+  NP5,XSHIFT,YSHIFT, NAIRFO, XI(202),YI(202), 

+  COSTHL(201),SINTHL(201) 
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110 


COMMON  /COF/  A(201,211),KEQNS 

COMMON  /SING/  Q( 200) ,GAMMA( 2) ,QK( 200) ,GAMK(2) 

COMMON  /WAK/  VYW(2) ,VXW(2) ,WAKE( 2) ,DT 

COMMON  /CORV/  CV(  2  ,  200)  ,XC(  2  ,200)  ,YC(2,200)  ,M,TD,CWX(2 ,200) , 
+  CVVY(2,200),XCI(2,200),YCI(2,200) 

COMMON  /INF1/  AAN(201,201),BBN(201,201),AYNP1(2,201),BYNF1(2,201) 
+  SUMAAN(201,2),SUMBBN(201,2) 

COMMON  /INF2/  SUMCCNC 201) ,SUMCCT( 201) ,CYNP1( 2,200) , 
+  CXNP1(2,200) 

COMMON  /GUST/  UG( 200) ,VG( 200) ,XGF,UGUST,VGUST 

COMMON  /MATRIX/  AMAT( 20 1,211) 

COMMON  /GEOM/  SINALF(2) ,COSALF(2) ,OMEGA(2) ,UX(2) ,UY(2) , PIVOT, 
+  XPRM,YPRM 

NEQS    =  NODTOT*  NAIRFO 

IF  (NITR  .GT.  0)  GO  TO  91 

SET/RESET  LHS  MATRIX  A(I,J)  FOR  EACH  NEW  TIME  STEP 

DO  110   I  =  l,2*NODTOT 
DO  110   J  =  l,2*NODTOT 
A(I,J)  =  AAN(I,J) 
IF  (M.NE.O)   GOTO  91 

FILL  IN  THE  RIGHT  HAND  SIDE  FOR  STEADY  FLOW 

DO  310  L=l, NAIRFO 


LI 

KI 

DO  310 

XMID 

YMID 


I 


(L-l)*NODTOT 
(L-1)*NP1 
=  l,NODTOT 


0.5  *  (XKI+KT)  +  XKI+KI+1)) 
0.5  *  (YKI+KT)  +  YKI+KI+1)) 
-SUMBBN(I+LI,1) 
-SUMBBN(I+LI,2) 

310  A(I+LI,NP5)  =  SINTHL(I+KI)*COSALF(L)  -  COSTHL( I+KI)*SINALF(L) 
RETURN 


A(I+LI,NP3)  = 
A(I+LI,NP4)  = 


FILL  IN  THE  RIGHT  HAND  SIDE  FOR  UNSTEADY  FLOW 


91 


210 


I  = 
=  0. 
=  0. 


L=l, NAIRFO 
=  (L-l)-NODTOT 
=  (L-1)*NP1 
l,NODTOT 

5   *   (XKI+KI)   +  XKI+KI+1)) 
5   *   (YI(I+KI)   +  YKI+KI+1)) 

=  -SUMBBN(I+LI,1)  +  BBN(I+LI,NP3)*SS(1)/WAKE(1) 
=  -SUMBBN(I+LI,2)  +  BBN( I+LI ,NP4)*SS(2)/WAKE(2) 
=  -BBN( I+LI ,NP3)*GAMMA( 1)*SS( 1)/WAKE( 1)-BBN( I+LI ,NP4) 
+*GAMMA(2)*SS(2)/WAKE(2)  +  SINTHL( I+KI)*  (( 1. +UG( I+LI) )*COSALF(L)- 
+VG(I+LI)*SINALF(L)+UX(L))  -  COSTHL( I+KI)*  (( 1. +UG( I+LI)) 
+*SINALF( L)+VG( I+LI )*COSALF( L)+UY( L) )+OMEGA( L)*( YMID*SINTHL( I+KI ) 
+  +  XMID*COSTHL(I+KI)) 

ADD  CORE  VORTEX  CONTRIBUTION 


DO 

LI 

KI 

DO  210 

XMID 

YMID 

A(I+LI,NP3) 

A(I+LI,NP4) 

A(I+LI,NP5) 


IF  (M  .EQ.  1)  GOTO  210 
A(I+LI,NP5)  =  A(I+LI,NP5) 


SUMCCNC I+LI) 
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210  CONTINUE 

DO  300  I=l,2*NODTOT 
DO  300  J=l,2*N0DT0T+3 
300  AMAT(I,J)  =  A(I,J) 
RETURN 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C     SUBROUTINE  KUTTA  ( ALPHA, S I NALF, COS ALF, OMEGA ,UX,UY)  C 

c  c 

C  USING  KUTTA  CONDITION  TO  DETERMINE  VORTICITY  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 
c 
c 
c 
c 
c 
c 


GAMMA 

GAMK 
QK 


c 
c 


CIRCULATION  FOR  THE  STEADY  FLOW  CASE  AND  ALSO 
CIRCULATION  FOR  THE  PREVIOUS  TIME  STEP  FOR  UNSTEADY  CASE 
CIRCULATION  AT  THE  CURRENT  TIME  STEP  FOR  UNSTEADY  CASE 
SOURCE  STRENGTH  AT  THE  CURRENT  TIME  STEP 
QK(I)  =  B1(I)*GAMK(1)+B2(I)*GAMK(2)+B3(I) 

B1(I)  ...  THAT  PART  OF  THE  SOURCE  STRNGTH  WHICH  IS  INDUCED  BY  A 
CIRCULATION  OF  UNIT  STRENGTH  FROM  AIRFOIL  1 

B2(I)  ...  THAT  PART  OF  THE  SOURCE  STRNGTH  WHICH  IS  INDUCED  BY  A 
CIRCULATION  OF  UNIT  STRENGTH  FROM  AIRFOIL  2 

B3(I)  ...  THAT  PART  OF  THE  SOURCE  STRNGTH  WHICH  IS  INDEPENDENT 
OF  THE  CIRCULATION  FROM  BOTH  AIRFOILS 

AA(I,J)  .  THAT  PART  OF  THE  TANGENTIAL  VELOCITY  AT  THE  TRAILING 

EDGE  PANELS  WHICH  IS  INDUCED  BY  A  CIRCULATION  OF  UNIT 
STRENGTH  BY  AIRFOIL  J 

BB(I)  ...  THAT  PART  OF  THE  TANGENTIAL  VELOCITY  AT  THE  TRAILING 
EDGE  PANELS  WHICH  IS  INDEPENDENT  OF  THE  CIRCULATION 

SUBROUTINE  KUTTA( NITR , PVTAG) 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y( 202) , 
+  C0STHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 

+  NP5 ,XSHIFT,YSHIFT,NAIRFO,XI(202) ,YI(202) , 

+  COSTHL(201),SINTHL(201) 

COMMON  /COF/  A(201,211),KEQNS 

COMMON  /SING/  Q( 200) ,GAMMA( 2) ,QK( 200) ,GAMK( 2) 

COMMON  /WAK/  VYW( 2) , VXW( 2) ,WAKE( 2) ,DT 

COMMON  /CORV/  CV( 2 , 200) ,XC( 2 , 200) , YC( 2 ,200) ,M,TD,CVVX( 2 ,200) , 
+  CVVY(2,200),XCI(2,200),YCI(2,200) 

COMMON  /INF1/  AAN(201,201),BBN(201,201),AYNP1(2,201),BYNP1(2,201), 
+  SUMAAN(201,2),SUMBBN(201,2) 

COMMON  /INF2/  SUMCCN(201) ,SUMCCT(201) ,CYNP1(2,200) , 
+  CXNP1(2,200) 

COMMON  /GUST/  UG( 200) ,VG( 200) ,XGF,UGUST,VGUST 

COMMON  /MATRIX/  AMAT( 201 ,211) 

COMMON  /BMAT/  B1,B2,B3 

COMMON  /GEOM/  SINALF( 2) ,COSALF(2) ,OMEGA( 2) ,UX(2) ,UY(2) .PIVOT, 
+  XPRM.YPRM 

COMMON  /GIES/  NGIES 

DIMENSION  BK200)  ,B2(200)  ,B3(200)  ,AA(4,2)  ,BB(5)  ,EE(2)  ,FF(2), 
+GG(2),RADI(2),AAA(2),BBB(2),CCC'(2),DDD(2),EEE(2),FFF(2), 
+DG(2),E(2),GAMA(2),DELG(2),DGAMK(2,3),COEFL(2,3) 

RETRIEVE  SOLUTION  FROM  A -MATRIX  FOR  THE  SOURCE  STRENGTH  IN  TERMS 
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50 


419 
430 
425 


THE  TWO  GAMK  CONTRIBUTION  AND  THE  CONSTANT  COMPONENT 

COUNT   =  0 

DO  50    I  =  l,NODTOT*2 
B1(I)   =  A(I,NP3) 
B2(I)   =  A(I,NP4) 
B3(I)    =  A(I,NP5) 
IF  (M.  GT.  0)   GOTO  400 


450 


460 


470 


STEADY  KUTTA  CONDITION 


430 
(K  . 
(K  . 


K 
EQ. 
EQ. 


COMPUTE 

DO  425 

LI 

KI 

KK 

DO 

IF 

IF 

XMID    =  0 

YMID    =  0 

AA(K+KK,1) 

AA(K+KK,2) 

BB(K+KK)  = 

DO  419   J  = 

AA(K+KK,1) 

AA(K+KK,2) 

BB(K+KK) 

CONTINUE 

CONTINUE 

CONTINUE 


TANGENTIAL  VELOCITIES 
L  =  l,NAIRFO 
=  (L-l)*NODTOT 
=  (L-1)*NP1 
(L-l)*2 
=  1,2 
1)  I  =  1 

I  =  NODTOT 
(XI(I+KI) 
(YI(I+KI) 
=  SUMAAN(I+LI 
=  SUMAAN(I+LI,2) 

COSALF(L)*COSTHL(I+KI)  +  SINALF(L)*SINTHL( I+KI) 
1 , 2*NODTOT 

=  AA(K+KK,1)  "  BBN(I+LI,J)*B1(J) 
=  AA(K+KK,2)  -  BBN(I+LI,J)*B2(J) 
=  BB(K+KK)  -  BBN(I+LI,J)*B3(J) 


2) 
.  5 

.  5 


+  XI(I+KI+1)) 
+  YI(I+KI+1)) 
,1) 


SET  UP  KUTTA  CONDITIONS 

DO  450  1=1,2 

LI      =  (I-l)*NAIRFO+l 
DGAMK(I,1)  =  AA(LI,1)  +  AA(LI+1,1) 
DGAMK(I,2)  =  AA(LI,2)  +  AA(LI+1,2) 
DGAMK(I,3)  =  -(BB(LI)  +  BB(LI+1)) 

SOLVE  KUTTA  CONDITION  BY  GAUSS 

R=DGAMK( 2,1) /DGAMK( 1,1) 

DO  460  K  =  2,3 

DGAMK(2,K)  =  DGAMK(2,K) -R*DGAMK( 1,K) 

GAMMA(2)  =  DGAMK(2,3)/DGAMK(2,2) 

GAMMA(l)  =  (DGAMK(1,3)-DGAMK(1,2)*GAMMA(2))/DGAMK(1,1) 

CALCULATE  SOURCE  STRENGTH 

DO  470  L  =  l,NAIRFO 

LI      =(L-l)*NODTOT 

DO  470   I  =  1, NODTOT 

Q(I+LI)   =  GAMMA(1)*B1(I+LI)  +  GAMMA(2)*B2( I+LI)  +  B3(I+LI) 

RETURN 
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c 

c 
c 

c 
c 


400 


119 


120 


c 

c 
c 


UNSTEADY  KUTTA  CONDITION 

FIND  VT  AT  TRAILING  EDGE  PANELS 


INITIAL  GUESS 
IF  (NITR 
GAMK(l)  : 
GAMK(2) 
END  IF 


FOR  GAMK 
EQ.  0)  THEN 
GAMMA(l) 
GAMMA(2) 


FROM  PREVIOUS  TIME  STEP 


DO 

LI 

KI 

KK 

DO 

IF 

IF 

XMID 

YMID 

VTANG 


125   L  =  l.NAIRFO 

=  (L-l)*NODTOT 
=  (L-1)*NP1 
=  (L-l)*2 
K  =  1,2 
EQ.  1)  I  =  1 
EQ.  2)  I  =  NODTOT 
=  0.5  *  (XI(I+KI) 
=  0.5  *  (YI(I+KI) 


130 
(K  . 
(K  . 


+  XKI+KI+1)) 
+  YI(I+KI+1)) 
( ( 1.  +UG( I+LI ) )*COSALF( L) -VG( I+LI )*SINALF( L)+UX( L) ) 
+*COSTHL(  I+KI )+  ( ( 1.  +UG(  I+LI )  )*SINALF(  L)+VG(  I+LI )*COSALF( L)+UY(  L)  ) 
+*SINTHL(I+KI)  +  OMEGA(L)*(YMID*COSTHL(I+KI) 
+-  XMID*SINTHL(I+KI)) 

AA(K+KK,1)   =  -  AAN(I+LI,NP3)*SS(1)/WAKE(1) 

AA(K+KK,2)   =  -  AAN(I+LI,NP4)*SS(2)/WAKE(2) 

BB(K+KK)   =  VTANG  +  AAN( I+LI ,NP3)*SS( 1)*GAMMA( 1)/WAKE( 1)  + 
+  AAN(I+LI,NP4)*SS(2)*GAMMA(2)/WAKE(2) 

DO  119   J  =  1, NODTOT 

AA(K+KK,1)   =  AA(K+KK,1)  +  AAN(I+LI,J)  -  BBN( I+LI , J)*B1( J) 

AA(K+KK,2)   =  AA(K+KK,2)  +  AAN( I+LI , J+NODTOT)  -  BBN( I+LI , J)*B2( J) 

BB(K+KK)   =  BB(K+KK)  -  BBN( I+LI , J)*B3( J) 

CONTINUE 

DO  120  J J  =  NODTOT+l,2*NODTOT 

AA(K+KK,1)   =  AA(K+KK,1)  -  BBN( I+LI , JJ)*B1( JJ) 

AA(K+KK,2)   =  AA(K+KK,2)  -  BBN( I+LI , JJ)*B2( JJ) 

BBCK+KK)   =  BB(K+KK)  -  BBN(  I+LI ,  JJ)*B3(  JJ) 

CONTINUE 

ADD  CORE  VORTEX  CONTRIBUTION 


100 
130 
125 


IF  (M.  LE. 

BB(K+KK) 

CONTINUE 

CONTINUE 

CONTINUE 

IF  (NGIES 


1)  GOTO  100 
=  BB(K+KK) 


+  SUMCCTCI+LI) 


EQ.  1)  GOTO  145 


SATISFYING  KUTTA  CONDITION 


SOLVE  FOR  VORTEX  STRENGTH 


DO  135  I  =  1,2 

LI  =  (I-l)*NAIRFO+l 

AAA(I)  =  AA(LI,1)**2-AA(LI+1,1)**2 

BBBCI)  =  AA(LI,2)**2-AA(LI+1,2)**2 

CCC(I)  =  2*(AA(LI,l)*BB(LI)-AA(LI+l,l)*BB(LI+l)-(2-I)*SS(l)/DT) 

DDD(I)  =  2*(AA(LI,2)*BB(LI)-AA(LI+1,2)*BB(LI+1)-(I-1)*SS(2)/DT) 

EEE(I)  =  2*(AA(LI,1)*AA(LI,2)-AA(LI+1,1)*AA(LI+1,2)) 

FFF( I )  =  BB( LI )**2-BB(LI+l)**2+2*SS( I ) *GAMMA( I ) /DT 
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135   CONTINUE 
60  DO  140  I  =  1,2 

DGAMK(I,1)  =  2*AAA(I)*GAMK(1)+EEE(I)*GAMK(2)+CCC(I) 
DGAMK(I,2)  =  2*BBB(I)*GAMK(2)+EEE(I)*GAMK(1)+DDD(I) 
140  DGAMK(I,3)  =  -( AAA( I)*GAMK( 1)**2+CCC( I)*GAMK( 1)+BBB( I)* 
+  GAMK(2)**2+DDD(I)*GAMK(2)+EEE(I)*GAMK(1)*GAMK(2)+FFF(I)) 
C140  WRITE  (6,*)'DGAMK' ,I,DGAMK(I,1),DGAMK(I,2),DGAMK(I,3) 
DO  144  1=1,2 
DO  144  KKK=1,3 

144  COEFL(I,KKK)  =  DGAMK(I,KKK) 
C 

C     GAUSSIAN  ELIMINATION 
C 

R=DGAMK( 2,1) /DGAMK( 1,1) 

DO  200  K  =  2,3 
200  DGAMK(2,K)  =  DGAMK( 2 ,K) -R*DGAMK( 1,K) 

DG(2)  =  DGAMK(2,3)/DGAMK(2,2) 

DG(1)  =  (DGAMK(1,3)-DGAMK(1,2)*DG(2))/DGAMK(1,1) 

GAMK(l)  =  GAMK(1)+DG(1) 

GAMK(2)  =  GAMK(2)+DG(2) 

IF  ((ABS(DG(1)/GAMK(1))  .  LT.  .0001)  .AND. 
+    (ABS(DG(2)/GAMK(2)).LT.  .0001))   GO  TO  300 

COUNT  =  COUNT  +  1 

IF  (COUNT  .EQ.  50)  GO  TO  310 

GO  TO  60 
310  WRTTE(6,*)'KUTTA  CONDITION  NOT  SATISFIED  --  COUNT  EXCEEDED  50' 

RETURN 
C 

C   SET  UP  KUTTA  CONDITIONS  FOR  GIESING  CONDITION 
C 

145  DO  451  1=1,2 

LI      =  (I-1)*NAIRF0+1 

DGAMK(I.l)  =  AA(LI,1)  +  AA(LI+1,1) 

DGAMK(I,2)  =  AA(LI,2)  +  AA(LI+1,2) 
451  DGAMK(I,3)  =  -(BB(LI)  +  BB(LI+1)) 
C 

C   SOLVE  KUTTA  CONDITION  BY  GAUSS 
C 

R=DGAMK( 2,1) /DGAMK( 1,1) 

DO  461  K  =  2,3 
461  DGAMK(2,K)  =  DGAMK(  2  ,K)  -R*DGAMK(  1,K) 

GAMK(2)  =  DGAMK(2,3)/DGAMK(2,2) 

GAMK(l)  =  (DGAMK(1,3)-DGAMK(1,2)*GAMK(2))/DGAMK(1,1) 
C 

C     CALCULATE  SOURCE  STRENGTH 
C 
300  DO  160   L  =  l,NAIRFO 

LI      =(L-l)*NODTOT 

DO  160  I  =  l,NODTOT 
160  QK(I+LI)   =  GAMK(1)*B1(I+LI)  +  GAMK(2)*B2(I+LI)  +  B3(I+LI) 
C 

C     CALCULATE  TANGENTIAL  VELOCITY  AT  THE  TRAILING  EDGE  BY  BACK 
C     SUBSTITUTION.  THE  PRODUCT  OF  THE  TRAILING  EDGE  VELOCITY 
C     SHOULD  BE  NEGATIVE  .  THIS  IS  CHECKED  IN  THE  MAIN  PROGRAM. 
C 

VTAG1   =  AA(1,1)*GAMK(1)+AA(1,2)*GAMK(2)+BB(1) 


120 


VTAG50  =  AA(2,1)*GAMK(1)+AA(2,2)*GAMK(2)+BB(2) 
VTAG51  =  AA(3,1)*GAMK(1)+AA(3,2)*GAMK(2)+BB(3) 
VTA100  =  AA(4,1)*GAMK(1)+AA(4,2)*GAMK(2)+BB(4) 
PVTAG   =  VTAG1*VTAG50 
2222  F0RMAT(1X,'VTAG' ,4F10.  6) 
RETURN 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C     SUBROUTINE  TEWAK  (SINALF,COSALF)  C 

c  c 

C  COMPUTE  WAKE  ELEMENT  AT  THE  TRAILING  EDGE  C 

C  C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C     VYVK(L)  Y-VELO  COMP  AT  THE  MID-POINT  OF  THE  WAKE  ELEMENT 

C  WITH  RESPECT  TO  THE  LOCAL  FROZEN  FRAME  OF  REFERENCE 

C     VXWK(L)  X-VELO  COMP  AT  THE  MID-POINT  OF  THE  WAKE  ELEMENT 

C                 WITH  RESPECT  TO  THE  LOCAL  FROZEN  FRAME  OF  REFERENCE 
C 

SUBROUTINE  TEWAK 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X( 202) , Y( 202) , 
+  COSTHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 

+  NP5,XSHIFT,YSHIFT,NAIRFO,XI(202),YI(202), 

+  COSTHL(201),SINTHL(201) 

COMMON  /COF/  A(201,211),KEQNS 

COMMON  /SING/  Q( 200) ,GAMMA( 2) ,QK( 200) ,GAMK( 2) 

COMMON  /WAK/  VYW( 2) , VXW( 2) ,WAKE( 2) ,DT 

COMMON  /WAK 2/  VYWK( 2) ,VXWK(2) 

COMMON  /CORV/  CV( 2 ,200) ,XC( 2,200) ,YC(2,200) ,M,TD,CVVX(2 ,200) , 
+  CVVY(2,200) ,XCI(2,200),YCI(2,200) 

COMMON  /INF1/  AAN(201,201),BBN(201,201),AYNP1(2,201),BYNP1(2,201), 
+  SUMAAN(201,2),SUMBBN(201,2) 

COMMON  /INF2/  SUMCCN( 201) ,SUMCCT(201) ,CYNP1( 2,200) , 
+  CXNP1(2,200) 

COMMON  /GEOM/  SINALF( 2) ,COSALF(2) ,OMEGA( 2) ,UX(2) ,UY( 2) , PIVOT, 
+  XPRM , YPRM 

COMMON  /GUST/  UG( 200) ,VG( 200) ,XGF,UGUST, VGUST 

DIMENSION  C0NT1X(2),C0NT1Y(2),C0NT2X(2),C0NT2Y(2), 
+  C0NT3X(2),C0NT3Y(2),C0NT4X(2),C0NT4Y(2),C0NT5X(2),C0NT5Y(2) 


10 


2221 


DO  20  I  = 

l,NAIRFO 

NN 

(I-l)-NPl 

XMID 

0.5  *  (X(NP1+NN)  + 

X(NP2+I-1)) 

YMID 

0.5  *  (Y(NP1+NN)  + 

Y(NP2+I-1)) 

UGW 

0.0 

VGW 

0.0 

XG 

XMID 

IF  (XG  .GT.  XGF)  GO  TO  10 

UGW 

UGUST 

VGW 

VGUST 

VYWK(I)  = 

(l.+UGW)*SINALF(I)+VGW*COSALF(I) 

VXWK(I)  = 

(l.+UGW)*COSALF(I)- 

-VGW*SINALF(I) 

FORMAT( IX 

, 'SINALF  : ' ,4F10.5) 

CONTIY(I) 

=  ( 1.  +UGW)*SINALF( I)+VGW*COSALF( I) 

CONTIX(I) 

=  ( 1.  +UGW)*COSALF(  I )  -VGW*SINALF(  I ) 

CONT2Y(I) 

=  0 
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C0NT2X(I)  =  0 
C0NT3YU)  =  0 
C0NT3X(I)  =  0 

CONTRIBUTION  FROM  THE  SOURCE  AND  WAKE  DIST  FROM  THE  AIRFOIL 


120 

110 


125 


DO  110  L  =  l.NAIRFO 

KJ      =  (L-l)*NODTOT 

DO  120  J  =  l,NODTOT 

VYWK(I)  =  VYWK(I)  +  AYNP1(I,J+KJ)*QK(J+KJ)  + 

VXVK(I)  =  VXWK(I)  -  BYNP1(I,J+KJ)*QK(J+KJ)  + 

CONT2Y(I)  =  CONT2Y(I)+  AYNP1( I , J+KJ)*QK( J+KJ) 

CONT2X(I)  =  CONT2X(I)-  BYNP1( I , J+KJ)*QK( J+KJ) 

CONT3Y(I)  =  CONT3Y(I)+  BYNP1( I , J+KJ)*GAMK(L) 

CONT3X(I)  =  CONT3X(I)+  AYNP1( I , J+KJ)*GAMK(L) 

CONTINUE 


ADD  WAKE  ELEMENT  CONTRIBUTION  OF  THE  OTHER  AIRFOIL 

J       =  3-1 

VYWK(I)  =  VYWK(I)+  BYNP1(I,NP3)*(GAMMA(J)-GAMK(J)) 
+  *  SS(J)/WAKE(J) 

VXWK(I)  =  VXWK(I)+  AYNP1(I,NP3)*(GAMMA(J)-GAMK(J)) 
+  *  SS(J)/WAKE(J) 

C0NT4Y(I)  =  BYNP1(I,NP3)*(GAMMA(J)-GAMK(J)) 
+  *  SS(J)/WAKE(J) 

C0NT4X(I)  =  AYNP1(I,NP3)*(GAMMA(J)-GAMK(J)) 
+  *  SS(J)/WAKE(J) 


BYNP1( I , J+KJ)*GAMK( L) 
A YNP 1 ( I , J+K J ) *G AMK ( L ) 


c 
c 
c 

ADD  CORE  VORTEX  CONTRIBUTION 

CONT5Y(I)  =  0 

CONT5X(I)  =  0 

IF  (M  .EQ.  1)  GO  TO  140 

MMl     =  M  -  1 

DO  150  L  =  1,NAIRF0 

KN      =(L-1)*MM1 

DO  130  N  =  1,MM1 

VYWK(I)  =  VYWK(I)  +  CYNP1(I,N+KN)*CV(L,N) 

VXWK(I)  =  VXWK(I)  +  CXNP1(I,N+KN)*CV(L,N) 

CONT5Y(I)  =  CYNP1(I,N+KN)*CV(L,N) 

130 

C0NT5X(I)  =  CXNP1(I,N+KN)*CV(L,N) 

150 

CONTINUE 

140 

CONTINUE 

20 

CONTINUE 

2222 

FORMAT  (IX,  'VXWK  :',6F10. 5) 

2223 

FORMAT  (IX,  'VYWK  :',6F10. 5) 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c                                            c 

C     SUBROUTINE  PRESS  (SINALF,COSALF, OMEGA, UX,UY)                   C 
C                                                               C 
C            COMPUTE  UNSTEADY  FLOW  PRESSURE  DISTRIBUTION             C 

C 

AND  VELOCITY  POTENTIAL  AT  MID-POINTS  OF  PANELS      C 

C 

C 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C  VTANG  TANGENTIAL  VELO  COMP  OF  A  FIXED  POINT  OF  THE  MOVING 

C  LOCAL  FRAME  OF  REFERENCE. 

PHIK(I) DISTURBANCE  VELO  POTENTIAL  AT  THE  MID  POINT  OF  THE 

C  I-TH  PANE  AT  THE  CURRENT  TIME  STEP  FOR  THE  UNSTEADY 

C  FLOW  CASE. 

C  PHI(I)  DISTURBANCE  VELO  POTENTIAL  AT  THE  MID  POINT  OF  THE 

C  I-TH  PANEL  FOR  THE  PREVIOUS  TIME  STEP 

C  PHILE(L) DIFFERENCE  OF  THE  POTENTIALS  OF  THE  LEADING  EDGE  TO 

C  THE  LOWER  TRAILING  EDGE  FOR  THE  RESPECTIVE  AIRFOIL 

C  PINK(L)  DIFFERENCE  OF  THE  POTENTIALS  AT  A  POINT  100  CHORD 

C  LENGTH  UPSTREAM  OF  THE  LE  FOR  THE  RESPECTIVE  AIRFOIL 

C  SUMC(L)  GAMMA  ASSOCIATED  WITH  THE  INTEGRATION  OF  THE  DISTURB- 

C  ANCE  VELOCITY  AROUND  THE  WHOLE  AIRFOIL 

C 

C  THE  FOLLOWING  INFLUENCE  COEF  ARE  COMPUTED  FOR  A  FINER  GRID  ON  THE 

C  AIRFOIL  SO  AS  TO  OBTAIN  A  MORE  ACCURATE  VELO  POTENTIAL  AT  THE  TE 

C  THE  INFLUENCE  COEF  ON  THE  I-TH  PANEL  FROM  THE  J-TH  PANEL  OF  THE 

C  AIRFOIL  REMAINS  THE  SAME  FOR  ALL  TIME  STEP 

C  

c 

C  AANP1(I,J,K). .  NORMAL  VELOCITY  INDUCED  AT  THE  I-TH  PANEL  SUB  NODE  K 

C  OF  AIRFOIL  1  DUE  TO  UNIT  STRENGTH  DIST  SOURCE  STRENG- 

C  TH  ON  THE  J-TH  PANEL  OF  AIRFOIL  1 

C  BBNP1(I,J,K)..  NORMAL  VELOCITY  INDUCED  AT  THE  I-TH  PANEL  SUB  NODE  K 

C  OF  AIRFOIL  1  DUE  TO  UNIT  STRENGTH  DIST  VORTICITY  STR- 

C  ENGTH  ON  THE  J-TH  PANEL  OF  AIRFOIL  1 

C  AANP2(I,J,K). .  NORMAL  VELOCITY  INDUCED  AT  THE  I-TH  PANEL  SUB  NODE  K 

C  OF  AIRFOIL  2  DUE  TO  UNIT  STRENGTH  DIST  SOURCE  STRENG- 

C  TH  ON  THE  J-TH  PANEL  OF  AIRFOIL  2 

C  BBNP2(I,J,K). .  NORMAL  VELOCITY  INDUCED  AT  THE  I-TH  PANEL  SUB  NODE  K 

C  OF  AIRFOIL  2  DUE  TO  UNIT  STRENGTH  DIST  VORTICITY  STR- 

C  ENGTH  ON  THE  J-TH  PANEL  OF  AIRFOIL  2 

C 

C 

c 

c    

SUBROUTINE  PRESS 

COMMON  /BOD/  IFLAG ,NLOWER,NUPPER ,NODTOT,X(  202) , Y(  202)  , 
+  COSTHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 

+  NP5,XSHIFT,YSHIFT,NAIRFO,XI(202),YI(202), 

+  COSTHL(201),SINTHL(201) 

COMMON  /CPD/  CP(200),SCL,T,SCM,SGAM 

COMMON  /NUM/  PI,PI2INV 

COMMON  /SING/  Q( 200) ,GAMMA(2) ,QK(200) ,GAMK(2) 

COMMON  /WAK/  VYW(2) , VXW( 2) ,WAKE(2) ,DT 

COMMON  /CORV/  CV(  2  ,  200)  ,XC(  2  ,200)  ,  YC(  2 ,200)  ,M,TD,CWX(2 ,200) , 
+  CVVY(2,200),XCI(2,200),YCI(2,200) 

COMMON  /INF1/  AAN(201,201),BBN(201,201),AYNP1(2,201),BYNP1(2,201), 
+  SUMAAN(201,2),SUMBBN(201,2) 

COMMON  /INF2/  SUMCCN( 201) ,SUMCCT( 201) ,CYNP1( 2,200) , 
+  CXNP1(2,200) 

COMMON  /INFL3/  AANP( 101 , 101,6) ,BBNP( 101 , 101,6) 

COMMON  /POT/  PHI(200),PHIK(200) 

COMMON  /GUST/   UG(200) ,VG(200) ,XGF,UGUST,VGUST 
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COMMON  /EXTV/  UE( 200) ,VN(200) 

COMMON  /BMAT/  B1,B2,B3 

COMMON  /GEOM/  SINALF( 2) ,COSALF( 2) ,0MEGA(2) ,UX(2) ,UY(2) , PIVOT, 
+  XPRM.YPRM 

DIMENSION  PHITEL(2),PHITEU(2),WGHT(5),PLOC(5),SUMC(2), 
+AANP1(50,50,9),AANP2(50,50,9),BBNP1(50,50,9),BBNP2(50,50,9), 
+COSTHP( 102,6) ,SINTHP( 102,6) ,AANP4(2) ,PHIL( 102) ,UGU( 100,6) , 
+VGU(100,6),PDUM(2),DP(2),BBNP4(2),B1(200))B2(200),B3(200),PHILE(2) 
3333  DIMENSION  CON2(2) ,VTGT(2,50) ,WAKCON(2,50) ,VOTCON(2,50) , 

+        DUMCON(2,50) 
7733  DIMENSION  CP1(200) 
4444  REAL  *8  PINK(2),VELX 

DATA  WGHT/.  11846344,. 
+  .23931434,. 

DATA  PLOC/.  04691008,. 
+  .76923466,. 

WRITE  (6,1000) 

IF  (M  .GT.  1)   GO  TO 


23931434,. 
11846344/ 
23076535,. 
95300899/ 

510 


28444444, 
50000000, 


100 


COMPUTE  THE  INFLUENCE  COEFFICIENT  FOR  A  FINER  GRID 

ON  THE  AIRFOIL  BY  THE  SAME  AIRFOIL  .  .  . 

(  ONLY  COMPUTED  ONCE  ) 


200 
200 


DO  200  L 

LI 

KI 

DO 

DO 

DX 

DY 

XMID 

YMID 

DO  200 

FLOG 

FTAN 

IF  (J+LI 

DXJ 

DXJP 

DYJ  ' 

DYJP 

FLOG 

FTAN 

CTIMTJ 


=  l,NAIRFO 

(L-l)*NODTOT 

(L-l)--NPl 

=  l,NODTOT 

=  1,5 

PLOC(K)*(X(I+KI+l)-X(I+KI)) 

PLOC(K)*(Y(I+KI+l)-Y(I+KI)) 

X(I+KI)  +  DX 

Y(I+KI)  +  DY 

=  l,NODTOT 

0.0 

PI 

EQ.  I+LI)     GO  TO  100 

XMID  -  X(J+KI) 

XMID  -  X(J+KI+1) 

YMID  -  Y(J+KI) 

YMID  -  Y(J+KI+1) 

.  5*ALOG( (DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ) ) 

ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 


COSTHE(I+KI)*COSTHE(J+KI) 

SINTHE(I+KI)*SINTHE(J+KI) 

STIMTJ  =  SINTHE(I+KI)*COSTHE(J+KI) 

COSTHE(I+KI)*SINTHE(J+KI) 

2)  GOTO  120 

PI2INV*(FTAN*CTIMTJ  • 
PI2INV*(FLOG'>CTIMTJ 


120 

200 

510 


IF  (L  .EQ. 

AANP1(I,J,K) 

BBNP1(I,J,K) 

GOTO  200 

AANP2(I,J,K) 

BBNP2(I,J,K) 

CONTINUE 

CONTINUE 


PI2INV*(FTAN*CTIMTJ 
PI2INV*(FLOG*CTIMTJ 


FLOG*STIMTJ) 
FTAN*STIMTJ) 

FLOG^STIMTJ) 
FTAN*STIMTJ) 
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FIND  TANGENTIAL  VELOCITY  AT  THE  MIDPOINT  OF  THE  I-TH  PANEL 

DO  600  L  =  l,NAIRFO 
LI=(L-l)*NODTOT 
KI=(L-1)*NP1 
SL*MC(L)   =0.0 
DO  600   I  =  l,NODTOT 
CONTR1      =  0.  0 
CONTR2      =  0.  0 
WAKCON(L,I)  =  0.  0 
VOTCON(L,I)  =  0.  0 
DUMCON(L,I)  =  0.  0 


XMID  =  0.5  *  (XKI+KI)  +  XI(I+KI+1)) 

YMID  =  0.5  *  (YKI+KI)  +  YKI+KI+1)) 

DX  =  (XKI+KI+1)  -  XI(I+KI)) 

DY  =  (YKI+KI+1)  -  YKI+KI)) 

DIST  =  SQRT(DX*DX+DY*DY) 

ACCOUNT  FOR  THE  FREESTREAM  AND  THE  MOTION  OF  THE  AIRFOIL 

VSX     =  (l.+UG( I+LI ))*COSALF(L)-VG( I+LI )*SINALF(L)  +  OMEGA(L)* 
+  YMID  +  UX(L) 

VSY     =  (l.+UG( I+LI) )*SINALF(L)+VG( I+LI )*COSALF(L) 
+  -  OMEGA(L)*XMID  +  UY(L) 

VS      =  VSX*VSX  +  VSY*VSY 

VNORM   =  -VSX*SINTHL(I+KI)+VSY*COSTHL(I+KI) 

VTANG   =  ( ( ( 1 . +UG( I+LI ) )*COSALF( L) -VG( I+LI )*SINALF( L)+UX( L) ) 
+  *COSTHL(I+Kn+  ((l.+UG(  I+LI  ) )*SINALF(L)+VG(  I+LI )*COSALF(L)+UY(L) ) 
+  *SINTHL(I+KI)  +  OMEGA(L)*(YMID*COSTHL(I+KI) 
+   -  XMID*SINTHL(I+KI))) 

VTFREE   =  VTANG 

VACT    =  VTANG 

INTRODUCE  SMALLER  GRIDS  FOR  THE  PURPOSE  OF  THE  VELO  POTENTIAL. 
VELO  ONLY  CALCULATED  AT  THE  MIDPT  OF  THE  PANEL  WHERE  K  =  3 

DO  260  K  =  1,5 

DX  =  PLOC(K)*(X(I+KI+l)-X(I+KI)) 

DY  =  PLOC(K)*(Y(I+KI+l)-Y(I+KI)) 

XINT  =  X(I+KI)  +  DX 

TINT  =  YCI+KI)  +  DY 

VDUM  =0.0 

INFLUENCE  COEF  ON  I-TH  PANEL  DUE  TO  THE  WAKE  ELEMENTS 

DO  245  MM  =  1,  NAIRFO 

J  =  NP1*MM 

DXJ  =  XINT  -  X(J) 

DXJP  =  XINT  -  X(NAIRF0*NP1+MM) 

DYJ  =  YINT  -  Y(J) 

DYJP  =  YINT  -  Y(NAIRF0*NP1+MM) 

FLOG  =  . 5*ALOG((DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 

FTAN  =  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 

CTIMTJ  =  COSTHE(I+KI)^COSTHE(J)  + 
+  SINTHE(I+KI)*SINTHE(J) 

STIMTJ  =  SINTHE(I+KI)*COSTHE(J)  - 
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+  COSTHE(I+KI)*SINTHE(J) 

AANP4(MM)  =  PI2INV*(FTAN*CTIMTJ  +  FLOG*STIMTJ) 
245   BBNP4(MM)  =  PI2INV*(FL0G*CTIMTJ  -  FTAN*STIMTJ) 
C0NTR1  =  SS(1)*(GAMMA(1)-GAMK(1))*AANP4(1) 
+/WAKE(1)+SS(2)*(GAMMA(2)-GAMK(2))*AANP4(2)/WAKE(2) 
C 

C     CONTRIBUTION  TO  VELO  COMPONENT  BY  WAKE  ELEMENT. 
C 

IF  (K  .EQ.  3)  THEN 

VACT  =  VACT  +  SS(1)*(GAMMA(1)-GAMK(1))*AAN(I+LI, 
+NP3)/WAKE(1)+SS(2)*(GAMMA(2)-GAMK(2))*AAN(I+LI,NP4)/WAKE(2) 

VNORM  =  VNORM  +  SS(  1)*(GAMMA(  1)-GAMK(  1))*BBNP4(  1) 
+/WAKE(1)+SS(2)*(GAMMA(2)-GAMK(2))*BBNP4(2)/WAKE(2) 
END  IF 
C 
C  EFFECTS  ON  AIRFOIL  BY  THE  SAME  AIRFOIL 

C         

c 

C     INTEGRATION  AROUND  FIRST  AIRFOIL 

C     CONTRIBUTION  TO  VELO  COMP  BY  AIRFOIL  1  WHEN  K  =  3 

C 

DO  300  J  =  l.NODTOT 
IF  (L  .EQ.  2)  GOTO  270 

VDUM     =  VDUM-BBNP1(I,J,K)*QK(J)+AANP1(I,J,K)*GAMK(1) 
IF  (K  .EQ.  3)  THEN 
VACT  =  VACT-BBN(I+LI,J)*QK(J)+ 
+AAN(I+LI,J)*GAMK(1) 

VNORM  =  VN0RM+(AANP1(I,J,K)*QK(J))+ 
+(BBNP1(I,J,K)*GAMK(1)) 
END  IF 
GOTO  300 
C 

C     INTEGRATION  AROUND  SECOND  AIRFOIL 
C     CONTRIBUTION  TO  VELO  COMP  BY  AIRFOIL  2  WHEN  K  =  3 
C 
270  VDUM     =  VDUM-BBNP2(I,J,K)*QK(J+NODTOT)+AANP2(I,J,K)*GAMK(2) 
IF  (K  .EQ.  3)  THEN 

VACT  =  VACT-BBN(I+LI,J+NODTOT)*QK(J+NODTOT)+ 
+  AAN( I+LI , J+NCDTOT)*GAMK( 2 ) 

VNORM  =  VN0RM+(AANP2(I,J,K)*QK(J+N0DT0T))+ 
+  (BBNP2(I,J,K)*GAMK(2)) 
ENDIF 
300  CONTINUE 
C 
C  EFFECTS  ON  AIRFOIL  BY  THE  OTHER  AIRFOIL 

C        

c 

C     INTEGRATION  AROUND  BOTH  AIRFOIL 

C     CONTRIBUTION  TO  VELO  COMP  BY  BOTH  AIRFOIL  WHEN  K  =  3 

C 

DO  350  J  =  NODTOT+2,2*NODTOT+l 

DXJ     =  XINT  -  X(J-KI) 

DXJP    =  XINT  -  X(J-KI+1) 

DYJ     =  YINT  -  Y(J-KI) 

DYJP    =  YINT  -  Y(J-KI+1) 

FLOG    =  . 5*ALOG((DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 
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350 


FT AN    =  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 

CTIMTJ  =  COSTHE(I+KI)*COSTHE(J-KI)  + 
+  SINTHE(I+KI)*SINTHE(J-KI) 

STIMTJ  =  SINTHE(I+KI)*COSTHE(J-KI)  - 
+  COSTHE(I+KI)*SINTHE(J-KI) 

AANP3   =  PI2INV*(FTAN*CTIMTJ  +  FLOG*STIMTJ) 

BBNP3   =  PI2INV*(FLOG*CTIMTJ  -  FTAN*STIMTJ) 

IF  (K.  EQ.  3)  THEN 

VACT  =  VACT-BBN(I+LI,J-LI-1)*QK(J-LI-1)+ 
+        AAN(I+LI,J-LI-1)*GAMK(3-L) 

VNORM  =  VNORM+(AANP3*QK(J-LI-l))+ 
+        (BBNP3*GAMK(3-L)) 

ENDIF 

VDUM    =  VDUM  -  BBNP3*QK(J-LI-1)  +  AANP3*GAMK(3-L) 


c 
c 

CONTRIBUTION  BY  CORE  VORTICES 

c 

1.  TO  VELO  POTENTIAL 

c 
c 

2.  TO  VELO  COMPONENT  ONLY  WHEN  K  =  3 

IF  (M. EQ. 1)  GOTO  150 

MM1 

=  M  -  1 

SUMCN 

=  0.0 

SUMCT 

=  0.0 

DO  400 

MM  =  1,  NAIRFO 

KN 

=  (MM-1)*MM1 

DO  400 

N  =  1,MM1 

DXC 

=  XINT  -  XC(MM,N) 

DYC 

=  YINT  -  YC(MM,N) 

DISTC 

=  SQRT(DXC*DXC+DYC*DYC) 

COSTHN 

=  DXC/DISTC 

SINTHN 

=  DYC /DISTC 

CTIMTN 

=  COSTHE(I+KI)*COSTHN  +  SINTHE( I+KI)*SINTHN 

STIMTN 

=  SINTHE(I+KI)*COSTHN  -  COSTHE(I+KI)*SINTHN 

CCN 

=  -CTIMTN/DISTC 

CCT 

=  -STIMTN/DISTC 

SUMCN 

=  SUMCN  +  CCN-CV(MM,N) 

400 

SUMCT 

=  SUMCT  +  CCT*CV(MM,N) 

IF  (K  . 

EQ.  3)  THEN 

VACT  = 

VACT+SUMCT*PI2INV 

VNORM  = 

=  VNORM+(PI2INV-SUMCN) 

ENDIF 

CONTR2 

=  PI2INV-SUMCT 

150  CONTINUE 

VTANG   =  VTANG  +  DBLE( (CONTRl+CONTR2+VDUM)*WGHT(K) ) 

SUMC(L)  =  SUMC(L)+VDUM*DIST*WGHT(K) 

WAKCON(L,I)  =  WAKCON(L,I)+CONTRl 

VOTCON(L,I)  =  VOTCON(L,I)+CONTR2 

DUMCON(L,I)  =  DUMCON(L,I)+VDUM 
260  CONTINUE 

PHIK(I+LI)  =  ( (VTANG) -VTFREE)*DI ST 

CP(I+LI)        =  VS    -    (VACT*VACT) 
7755  CPKI+LI)     =  CP(I+LI) 

UECI+LI)        =  VACT 

VN(I+LI)        =  VNORM 
600     CONTINUE 
6688  FORMAT(2I5,7E15.  7) 
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3335  FORMAT  (2I5.7F10. 6) 

COMPUTE  DISTURBANCE  POTENTIAL  AT  THE  LEADING  EDGE  BY  LINE 
INTEGRAL  OF  THE  VELOCITY  FIELD 
FROM  UPSTREAM  (AT  INFINITY)  TO  THE  LEADING  EDGE 
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24 


TRY 

PINK(l) 
PINK(2) 
DO  56 

YMID 

XLE 

XL 

XL 

NPHI 

PDUM(L) 

DO  30 

FRACT 

XLP 

IF  (I  . 

XLP 

DELX 

XMID 

XMID 

YMID 

XL 

VELX 


(L-1)*YSHIFT 
•1)*XSHIFT 


=  1 
=  0.0 
=  0.  0 

L  =  l,NAIRFO 
=  PIVOT  *  SINALF(L)  + 
=  -PIVOT*COSALF(L)+(L- 
=  XLE 
=  0.0 

=  10  *  NLOWER 
=  0.0 
I  =  1,NPHI 

=  FLOAT(I)/FLOAT(NPHI) 

=  -100.0  *TRY  *  (1.0  -  COS(0.5*PI*FRACT))+XLE 
EQ.  1)  XLP  =  -0.000197+XLE 
=  -10.0  *  (1.0  -  COS(0.5*PI*FRACT)) 
=  XL  -  XLP 
=  0.5*(XL+XLP) 
=  0.5*(XL+XLP)*COSALF(L) 
=  0.5*(XL+XLP)*SINALF(L) 
=  XLP 
=  UGUST 


ADD  CONTRIBUTION  OF  J-TH  PANEL 


40 


20 
(J 


DO 

LJ 

KJ 

DO 

IF 

DXJ 

DXJP 

DYJ 

DYJP 

FLOG 

FTAN 

CALMTJ 

SALMTJ 

CALMTJ 

SALMTJ 

APY 

BPY 

VELX 

GO  TO 

DXJ 

DXJP 

DYJ 

DYJP 

FLOG 

FTAN 

CALMTJ 


MM  =  l,NAIRFO 

=(MM-l)*NODTOT 

=(MM-1)*NP1 

J  =  1.NP1 
EQ.  NP1)  GO  TO  24 

=  XMID  -  X(J+KJ) 

=  XMID  -  X(J+KJ+1) 

=  YMID  -  Y(J+KJ) 

=  YMID  -  Y(J+KJ+1) 

=  .  5*ALOG((DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 

=  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 

=  -COSTHE(J+KJ) 

=  SINTHECJ+KJ) 

=  -COSALF(L)*COSTHE(J+KJ)  -  SINALF(L)*SINTHE( J+KJ) 

=  -SINALF(L)*COSTHE(J+KJ)  +  COSALF(L)*SINTHE( J+KJ) 

=  PI2INV*(FTAN*CALMTJ  +  FLOG*SALMTJ) 

=  PI2INV*(FL0G*CALMTJ  -  FTAN*SALMTJ) 

=  VELX   -  DPROD(BPY,QK(J+LJ))  +  DPROD(GAMK(MM) ,APY) 


20 


XMID  -  X(J+KJ) 

XMID  -  X(2*NP1+MM) 

YMID  -  Y(J+KJ) 

YMID  -  Y(2*NP1+MM) 

.  5»ALOG( ( DXJP*DXJP+DYJP*DYJP) /( DXJ*DXJ+DYJ*DYJ) ) 

ATAN2( DYJP*DXJ-DXJP*DYJ ,DXJP*DXJ+DYJP*DYJ) 

-COSALF(L)*COSTHE(J+KJ)  -  SINALF( L)-SINTHE( J+KJ) 
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SALMTJ  =  -SINALF(L)*COSTHE(J+KJ)  +  COSALF(L)*SINTHE( J+KJ) 

APY     =  PI2INV*(FTAN*CALMTJ  +  FLOG*SALMTJ) 

DUMMY   =  SS(MM)*(GAMMA(MM)-GAMK(MM))/WAKE(MM) 

VELX    =  VELX  +  DPROD( APY, DUMMY) 
20   CONTINUE 
40   CONTINUE 
C 

C     ADD  CORE  VORTEX  CONTRIBUTION 
C 

IF  (M  .EQ.  1)  GO  TO  50 

MM1     =  M  -  1 

DO  70   II  =  1,NAIRF0 

KN      =  (II-1)*MM1 

DO  60   N  =  1,MM1 

DX      =  XMID  -  XC(II,N) 

DY      =  YMID  -  YC(II,N) 

DIST    =  SQRT(DX*DX+DY*DY) 

COSTHN  =  DX/DIST 

SINTHN  =  DY/DIST 

SALMTN  =  -SINALF(L)*COSTHN  +  COSALF(L)*SINTHN 

CPT     =  -PI2INV*SALMTN/DIST 
60   VELX    =  VELX  +  DPROD(CPT,CV( II ,N) ) 
70   CONTINUE 
50   CONTINUE 

PDUM(L)  =  PDUM(L)  +  VELX  *  DBLE(DELX) 
7771  FORMAT  (215 ,6E14. 6.2F10. 6) 
30   CONTINUE 

DP(L)  =  PDUM(L)-PINK(L) 
56   CONTINUE 

PINK(l)  =  PDUM(l) 

PINK(2)  =  PDUM(2) 
C 

C     COMPUTATION  OF  THE  VELOCITY  POTENTIAL  FOR  MIDPOINT  OF  EACH  PANEL 
C 

DO  240  L   =  l,NAIRFO 

LI         =  (L-l)*NODTOT 

PHP        =  -PINK(L) 
C 

C     BEGIN  WITH  LOWER  SURFACE 
C 

DO  230  I   =  NLOWER,l,-l 

PHC  =  PHP-PHIKd+LI) 

PHIK(I+LI)  =  0.5*(PHP+PHC) 
230   PHP       =  PHC 

PHITEL(L)  =  PHC 
C 

C     RESET  FOR  UPPER  SURFACE 
C 

PHP       =  -PINK(L) 

DO  250  I   =  NLOWER+l,NODTOT 

PHC       =  PHP+PHIK(I+LI) 

PHIK(I+LI)  =  0.5*(PHP+PHC) 
250   PHP       =  PHC 

PHITEU(L)  =  PHC 
240  CONTINUE 
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COMPUTE  CP  AT  MID  POINT  OF  I-TH  PANEL 


DO  295   L  =  l,NAIRFO 

LI      =  (L-l)*NODTOT 

KI      =  (L-1)*NP1 

SUMC(L)  =  SUMC(L)/SS(L) 

DO   290      I  =   1,N0DT0T 

XIMID        =   . 5*(XI(I+KI)   +  XKI+KI+1)) 

XMID  =  .5*(X(I+KI)   +  X(I+KI+1)) 

YMID  =   . 5*(Y(I+KI)   +  Y(I+KI+1)) 

CP(I+LI)=  CP(I+LI)    -   2.*(PHIK(I+LI)-PHI(I+LI))/DT 

WRITE  (6,1050)  I+LI, XIMID, XMID, YMID, QK( I+LI ),GAMK(L),CP( I+LI), 
+  UE( I+LI ) , VN( I+LI ) , PHIK( I+LI ) , PHI ( I+LI ) , SUMC( L) 
290  CONTINUE 

WRITE  (6,235)  PINK(L) 
295   CONTINUE 

235   FORMAT  (IX,' PHI  AT  LEADING  EDGE  =' ,F10. 6,/) 

1000  F0RMAT(/,4X,'J'  ,4X,'XI(J)'  ,5X,'X(J)'  ,6X,fY(J)'  ,6X,'Q(J)f  ,6X, 
+' GAMMA* ,4X.fCP(J)' .7X,'V(J)' ,6X,'VN(J)' ,5X,'PHIK' ,6X, 
+'PHI' ,  6X,fINTGAMKf ,/) 
1050  FORMAT(I5,12F10. 5) 
1200  FORMAT( IX, 'LENGTH  OF  LEADING  EDGE  INTEGRATION  IN  CHORD  =',F10.6/) 

RETURN 

END 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C     SUBROUTINE  CORVOR  ( SINALF,COSALF)  C 

c  c 

C  COMPUTE  THE  LOCAL  VELOCITIES  OF  CORE  VORTICES  C 

C  C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C  CCVX(I,N)  X-VELO  OF  THE  I-TH  AIRFOIL,  N-TH  CORE  VORTEX  WITH 

C  RESPECT  TO  THE  CURRENT  FROZEN  FRAME  OF  REFERENCE 

C  CCVY(I,N)  Y-VELO  OF  THE  I-TH  AIRFOIL,  N-TH  CORE  VORTEX  WITH 

C  RESPECT  TO  THE  CURRENT  FROZEN  FRAME  OF  REFERENCE 

C  XC(I,N)    X-COORD  OF  THE  LOCATION  OF  THE  I-TH  AIRFOIL,  N-TH 

C  CORE  VORTEX  W.  R.  T.  THE  GLOBAL  FRAME  OF  REFERENCE. 

C  YC(I,N)    Y-COORD  OF  THE  LOCATION  OF  THE  I-TH  AIRFOIL,  N-TH 

C  CORE  VORTEX  W. R. T.  THE  GLOBAL  FRAME  OF  REFERENCE. 

C  GAMY(I)    GLOBAL  Y-VELO  AT  A  LOCATION  OF  THE  I-TH  AIRFOIL 

C  CORE  VORTEX  DUE  TO  A  SOURCE  DIST  OF  UNIT  STRENGTH 

C  ON  ONE  PANEL 

C  GBMY(I)    GLOBAL  Y-VELO  AT  A 

C  CORE  VORTEX  DUE  TO 

C  ON  ONE  PANEL 

C  AMY(I)    LOCAL  Y-VELO  AT  A  LOCATION  OF  THE  I-TH  AIRFOIL 

C  CORE  VORTEX  DUE  TO  A  SOURCE  DIST  OF  UNIT  STRENGTH 

C  ON  ONE  PANEL 

C  BMY(I)    LOCAL  Y-VELO  AT  A  LOCATION  OF  THE  I-TH  AIRFOIL 

C  CORE  VORTEX  DUE  TO  A  VORTEX  DIST  OF  UNIT  STRENGTH 

C  ON  ONE  PANEL 

C  SUMAMY(I)  LOCAL  Y-VELO  AT  A  LOCATION  OF  THE  I-TH  AIRFOIL 

C  CORE  VORTEX  DUE  TO  A  SOURCE  DIST  OF  UNIT  STRENGTH 

C  ON  ALL  PANELS 


LOCATION  OF  THE  I-TH  AIRFOIL 
A  VORTEX  DIST  OF  UNIT  STRENGTH 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


c 
c 
c 


SUMBMY(I) 

GCMX 
GCMY 
CMX 
CMY 


LOCAL  Y-VELO  AT  A 
CORE  VORTEX  DUE  TO 
ON  ALL  PANELS 
GLOBAL  X-VELO  AT  A 
TO  ANOTHER  POINT 
GLOBAL  Y-VELO  AT  A 
TO  ANOTHER  POINT 
LOCAL  X-VELO  AT  A 
TO  ANOTHER  POINT 
LOCAL  Y-VELO  AT  A 
TO  ANOTHER  POINT 


LOCATION  OF  THE  I-TH  AIRFOIL 
A  VORTEX  DIST  OF  UNIT  STRENGTH 


LOCATION  OF  A  CORE  VORTEX  DUE 
VORTEX  OF  UNIT  STRENGTH 

LOCATION  OF  A  CORE  VORTEX  DUE 
VORTEX  OF  UNIT  STRENGTH 
LOCATION  OF  A  CORE  VORTEX  DUE 
VORTEX  OF  UNIT  STRENGTH 
LOCATION  OF  A  CORE  VORTEX  DUE 
VORTEX  OF  UNIT  STRENGTH 


SUBROUTINE  CORVOR 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X( 202) ,Y( 202) , 
+  COSTHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 

+  NP5,XSHIFT,YSHIFT,NAIRFO,XI(202),YI(202), 

+  COSTHL(201),SINTHL(201) 

COMMON  /SING/  Q( 200) ,GAMMA( 2) ,QK( 200) ,GAMK(2) 

COMMON  /WAK/  VYV( 2) , VXW( 2) ,WAKE( 2) ,DT 

COMMON  /CORV/  CV( 2 ,200) ,XC( 2 , 200) , YC( 2 , 200) ,M,TD,CCVX( 2 , 200) , 
+  CCVY(2,200),XCI(2,200),YCI(2,200) 

COMMON  /POT/  PHI(200),PHIK(200) 

COMMON  /GUST/  UG( 200) , VG( 200) ,XGF,UGUST, VGUST 

COMMON  /NUM/  PI.PI2INV 

COMMON  /GEOM/  SINALF( 2) ,COSALF( 2) ,OMEGA( 2) ,UX( 2) ,UY( 2) ,PIVOT, 
+  XPRM , YPRM 

DIMENSION  AMY( 2) ,BMY( 2) , SUMAMY( 2) ,SUMBMY( 2) , 
+VX(2),VY(2),GAMY(2),GBMY(2) 

IF  (M.  EQ.  1)  GOTO  40 

MM1     =  M  -  1 

ACCOUNT  FOR  THE  FREE  STREAM  INCLUDING  GUST  EFFECTS 


=  0.  0 
=  0.  0 

I  =  l.NAIRFO 
=  (I-1)*MM1 

N  =  1.MM1 
=  XC(I,N)*COSALF(I)  +  YC(I,N)*SINALF(I) 
.GT.  XGF)  GO  TO  5 
=  UGUST 
=  VGUST 


UGC 

VGC 

DO  15 

KN 

DO  10 

XG 

IF  (XG 

UGC 

VGC 

CONTINUE 

VY(I)   =  (l.+UGC)*SINALF(I)+VGC*COSALF(I) 

VX(I)   =  (l.+UGC)*COSALF(I)-VGC*SINALF(I) 

XMID   =  XC(I,N) 

YMID   =  YC(I,N) 

CALCULATE  THE  INFLUENCE  COEFFICIENT  DUE  TO  THE  AIRFOILS 

DO  25   L  =  l.NAIRFO 
LJ      =(L-l)*NODTOT 
KJ      =(L-1)*NP1 
SUMAMY(I)  =0.0 
SUMBMY(I)  =0.0 
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11 

12 


20 


25 


30 
35 


DO  20   J  =  1,NP1 

DXJ    =  XMID  -  X(J+KJ) 

DYJ    =  YMID  -  Y(J+KJ) 

IF  (J  .EQ.  NP1)  GO  TO  11 

DXJP   =  XMID  -  X(J+KJ+1) 

DYJP   =  YMID  -  Y(J+KJ+1) 

GO  TO  12 

DXJP 

DYJP 

FLOG 

FT  AN 

GAMY(I)  = 

GBMY(I)  =  PI2INV*(FL0G*C0STHE(J+KJ)  +  FTAN*SINTHE( J+KJ) ) 

AMY(I)   =  GAMY(I)*COSALF(I)-GBMY(I)*SINALF(I) 

BMY(I)   =  GAMY(I)*SINALF(I)+GBMY(I)*COSALF(I) 

IF  (J.  EQ.  NP1)  GOTO  20 

SUMAMY(I)  =  SUMAMY(I)  +  AMY(I) 

SUMBMY(I)  +  BMY(I) 

VY(I)  +  AMY(I)*QK(J+LJ) 

VX(I)  -  BMY(I)*QK(J+LJ) 


XMID  -  X(2*NP1+L) 
YMID  -  Y(2*NP1+L) 

. 5*AL0G( (DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ) ) 
ATAN2( DYJP*DXJ-DXJP*DYJ ,  DXJP*DXJ+DYJP*DYJ) 
PI2INV*(FTAN*C0STHE(J+KJ)  -  FLOG*SINTHE( J+KJ) ) 


SUMBMY(I)  = 
VY(I) 
VX(I) 
CONTINUE 


ADD  THE  CONTRIBUTION  DUE  TO  THE  WAKE  ELEMENTS 

VY(I)  =  VY(I)  +  SUMBMY(I)*GAMK(L)  +  SS(L)*BMY( I)* 
+        (GAMMA(L)-GAMK(L))/WAKE(L) 

VX(I)  =  VX(I)  +  SUMAMY(I)*GAMK(L)  +  SS(L)*AMY( I)* 
+         (GAMMA(L)-GAMK(L))/WAKE(L) 

CONTINUE 

CALCULATE  INFLUENCE  COEFFICIENT  DUE  TO  THE  WAKE 

DO  35  L  =  1,NAIRF0 

KMC  =  (L-1)*MM1 

DO  30  MC  =  1,MM1 

IF  ((L.EQ.  I)  .AND.  (MC.EQ.N))  GOTO  30 

DX  =  XMID  -  XC(L,MC) 

DY  =  YMID  -  YC(L,MC) 

DIST2  =  DX*DX+DY*DY 

GCMY  '  =  -PI2INV*DX/DIST2 

GCMX  =  +PI2INV*DY/DIST2 

CMY  =  GCMY»COSALF(I)+GCMX"SINALF(I) 

CMX  =  -GCMY*SINALF(I)+GCMX*COSALF(I) 

ADD  CORE  VORTEX  CONTRIBUTION 


VY(I)   = 
VX(I)  = 

CONTINUE 
CONTINUE 


VY(I)  +  CMY*CV(L,MC) 
VX(I)  +  C>DC^CV(L,MC) 


10 


CONVECTION  VELOCITY  OF  CORE  VORTICES  AT  NEXT  TIME  STEP 

CCVX(I,N)  =  VX(I) 
CCVY(I,N)  =  VY(I) 
CONTINUE 
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C  =  0 

C  =  1 

C  =2 

C  =3 


15   CONTINUE 
40   CONTINUE 
RETURN 
END 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  C 

C     SUBROUTINE  NEWPOSN  C 

c  c 

C     TRANSFORM  THE  RESPECTIVE  FROZEN  LOCAL  COORDINATES  AND  PANEL     C 
C     ANGLES  TO  THE  GLOBAL  FRAMES  OF  REFERENCE  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C     XSHIFT  GLOBAL  X  LOCATION  OF  THE  2-ND  AIRFOIL  PIVOT  POSITION 

C     YSHIFT  GLOBAL  Y  LOCATION  OF  THE  2-ND  AIRFOIL  PIVOT  POSITION 

ITYPE   FLAG  FOR  DEALING  WITH  THE  VARIOUS  COMPONENTS 

FOR  TRANSFORMING  AIRFOIL  COORDINATES 

FOR  TRANSFORMING  WAKE  ELEMENTS 

FOR  TRANSFORMING  MOST  RECENT  CORE  VORTEX  SHED 

FOR  TRANSFORMING  ALL  PREVIOUS  CORE  VORTICES 

XPRM   X  GLOBAL  RELATIVE  MOVEMENT  OF  THE  PIVOTS  POSITION 

C     YPRM   Y  GLOBAL  RELATIVE  MOVEMENT  OF  THE  PIVOTS  POSITION 

C     

C 

SUBROUTINE  NEWPOS( ITYPE) 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y( 202) , 
+  COSTHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 

+  NP5, XSHIFT, YSHIFT, NAIRFO,XI(202),YI(202), 

+  COSTHL(201),SINTHL(201) 

COMMON  /CORV/  CV( 2, 200) ,XC( 2 ,200) ,YC(2 ,200) ,M,TD,CVVX( 2,200) , 
+  CVVY(2,200),XCI(2,200),YCI(2,200) 

COMMON  /GEOM/  SINALF( 2) ,COSALF( 2) ,OMEGA( 2) ,UX( 2) ,UY(2) , PIVOT, 
+  XPRM, YPRM 

C 

C     CHECK  FLAG  TO  ASCERTAIN  TYPE  OF  TRANSFORMATION 
C 

IF  (  ITYPE  .EQ.  1  )  GOTO  30 
IF  (  ITYPE  .EQ.  2  )  GOTO  40 
IF  (  ITYPE  .EQ.  3  )  GOTO  45 
C 

C     AIRFOIL  COORDINATES  AND  LOCAL  ANGLES  TRANSFORMATION 
C 

DO  10  I  =  1,NP1 

X(I)    =  (XI(I)-PIVOT)*COSALF(l)+YI(I)*SINALF(l) 
Y(I)    =  (XI(I)-PIVOT)*(-SINALF(l))+YI(I)*COSALF(l) 
10  CONTINUE 

DO  20  I  =  NP1+1,2*NP1 

X(I)    =  (XI(I)-PIVOT)*COSALF(2)+YI(I)*SINALF(2)+ 
+         XSHIFT  +  XPRM 

Y(I)    =  (XI(I)-PIVOT)*(-SINALF(2))+YI(I)*COSALF(2)+ 
+  YSHIFT  +  YPRM 

20  CONTINUE 
C 

C     TRANSFORM  LOCAL  ANGLES 
C 

DO  222  L  =1,2 
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222 


30 


223 


40 


45 


46 
50 


K       =  (L-l)*(NODTOT+l) 
DO  222  J  =  1,N0DT0T 
DX      =  X(J+1+K)  -  X(J+K) 
DY      =  Y(J+1+K)  -  Y(J+K) 
DIST    =  SQRT(DX*DX  +DY*DY) 
SINTHE(J+K)  =  DY/DIST 
COSTHE(J+K)  =  DX/DIST 
GOTO  50 

WAKE  ELEMENT  TRANSFORMATION 

I  =  NP2 

X(I)  =  (XI(I)-PIV0T)*C0SALF(1)+YI(I)*SINALF(1) 

Y(I)  =  (XI(I)-PIV0T)*(-SINALF(1))+YI(I)*C0SALF(1) 

I  =  NP2  +  1 

X(I)  =  (XI(I)-PIV0T)*C0SALF(2)+YI(I)*SINALF(2)+ 
h  XSHIFT  +  XPRM 

Y(I)  =  (XI(I)-PIV0T)*(-SINALF(2))+YI(I)*C0SALF(2)+ 

V  YSHIFT  +  YPRM 

WAKE  ELEMENT  ANGLES  TRANSFORMATION 

DO  223  L  =1,2 

K       =  (L-1)*(N0DT0T+1) 

J       =  NP1 

DX      =  X(NP2+L-1)  -  X(NP1*L) 

DY      =  Y(NP2+L-1)  -  Y(NP1*L) 

DIST    =  SQRT(DX*DX  +DY*DY) 

SINTHE(J+K)  =  DY/DIST 

C0STHE(J+K)  =  DX/DIST 

GOTO  50 

MOST  RECENT  CORE  VORTEX  SHED  TRANSFORMATION 

XC(1,M)  =  (XCI(1,M)-PIV0T)*C0SALF(1)+YCI(1,M)*SINALF(1) 
YC(1,M)  =  (XCI(1,M)-PIV0T)*(-SINALF(1))+YCI(1,M)*C0SALF(1) 

XC(2,M)  =  (XCI(2,M)-PIV0T)*C0SALF(2)+YCI(2,M)*SINALF(2)+ 

V  XSHIFT  +  XPRM 

YC(2,M)  =  (XCI(2,M)-PIV0T)*(-SINALF(2))+YCI(2,M)*C0SALF(2)+ 

V  YSHIFT  +  YPRM 
GOTO  50 

ALL  PREVIOUS  CORE  VORTICES  TRANFORMATION 

DO  46   I  =  l.M-1 

XC(1,I)  =  (XCI(1,I)-PIV0T)*C0SALF(1)+YCI(1,I)*SINALF(1) 
YC(1,I)  =  (XCI(1,I)-PIV0T)*(-SINALF(1))+YCI(1,I)*C0SALF(1) 
XC(2,I)  =  (XCI(2,I)-PIVOT)*COSALF(2)+YCI(2,I)*SINALF(2)+ 

V  XSHIFT  +  XPRM 

YC(2,I)  =  (XCI(2,I)-PIV0T)*(-SINALF(2))+YCI(2,I)*C0SALF(2)+ 
»-  YSHIFT  +  YPRM 

CONTINUE 
CONTINUE 
RETURN 
END 
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APPENDIX      B.  EXAMPLE  LNPUT  DATA  FOR  PROGRAM  USPOTF2 

5 
************************************ 

AIRFOIL  TYPE  :  8.4%  THICKNESS  SYMMETRICAL  MISES  AIRFOIL 
AIRFOIL  COORDINATES  ARE  USER  INPUT  DATA  ( IFLAG  =  1) 
NLOWER  =  25  ,  NUPPER  =  25 

****************************************************************** 


01   25 

25 

2  0.0 

2.0 

1.  000000 

0. 994858 

0. 980866 

0. 

958884 

0. 

929536 

0.893455 

0. 851308 

0.803815 

0.  751753 

0. 

695948 

0. 

637271 

0.576620 

0.514918 

0.453098 

0. 392084 

0. 

332794 

0. 

276105 

0.222865 

0. 173861 

0.  129819 

0.091393 

0. 

059146 

0. 

033560 

0.015010 

0. 003767 

0. 000000 

0. 003767 

0. 

015008 

0. 

033560 

0.059146 

0. 091393 

0. 129819 

0.  173861 

0. 

222865 

0. 

276105 

0.332791 

0. 392082 

0.453095 

0.514915 

0. 

576617 

0. 

637266 

0. 695946 

0.  751750 

0. 803815 

0. 851308 

0. 

893455 

0. 

929536 

0.958884 

0. 980866 

0.  994858 

1.  000000 

0. 000000 

-0. 000782 

-0. 002784 

-0. 

005721 

-0. 

009351 

-0.013459 

-0. 017837 

-0.022285 

-0.026618 

-0. 

030671 

-0. 

034289 

-0.037341 

-0. 039712 

-0. 041314 

-0. 042083 

-0. 

041979 

-0. 

040979 

-0.039096 

-0.036360 

-0. 032820 

-0.028555 

-0. 

023651 

-0. 

018220 

-0.012379 

-0. 006259 

0. 000000 

0. 006259 

0. 

012379 

0. 

018220 

0.023651 

0. 028555 

0. 032820 

0. 036360 

0. 

039096 

0. 

040979 

0.041979 

0.042083 

0.  041314 

0.  039712 

0. 

037341 

0. 

034289 

0.030671 

0.026618 

0. 022285 

0.  017837 

0. 

013459 

0. 

009351 

0.005721 

0. 002784 

0. 000782 

0. 000000 

1.000000 

0.994858 

0. 980866 

0. 

958884 

0. 

929536 

0. 893455 

0.851308 

0.  803815 

0.  751753 

0. 

695948 

0. 

637271 

0.576620 

0.514918 

0. 453098 

0. 392084 

0. 

332794 

0. 

276105 

0.222865 

0. 173861 

0. 129819 

0. 091393 

0. 

059146 

0. 

033560 

0.015010 

0. 003767 

0. 000000 

0. 003767 

0. 

015008 

0. 

033560 

0.059146 

0.  091393 

0. 129819 

0. 173861 

0. 

222865 

0. 

276105 

0.332791 

0. 392082 

0.453095 

0. 514915 

0. 

576617 

0. 

637266 

0.695946 

0.  751750 

0. 803815 

0. 851308 

0. 

893455 

0. 

929536 

0.958884 

0. 980866 

0. 994858 

1.  000000 

0. 000000 

-0. 000782 

-0. 002784 

-0. 

005721 

-0. 

009351 

-0.013459 

-0. 017837 

-0. 022285 

-0. 026618 

-0. 

030671 

-0. 

034289 

-0.037341 

-0.039712 

-0. 041314 

-0.042083 

-0. 

041979 

-0. 

040979 

-0.039096 

-0. 036360 

-0. 032820 

-0.028555 

-0. 

023651 

-0. 

018220 

-0.012379 

-0. 006259 

0. 000000 

0. 006259 

0. 

012379 

0. 

018220 

0.023651 

0.028555 

0.032820 

0. 036360 

0. 

039096 

0. 

040979 

0.041979 

0. 042083 

0. 041314 

0. 039712 

0. 

037341 

0. 

034289 

0.030671 

0. 026618 

0. 022285 

0.017837 

0. 

013459 

0. 

009351 

0.005721 

0.002784 

0. 000782 

0. 000000 

0.0 

0.0 

-45.80000 

0. 

00    0 

.0 

. 0000      0 

0.  0 

0.00 

0.  0 

0. 

00    0 

.0 

.000       0.0 

0.0 

5.0000 

0.065 

0.001 

0. 

,00    4 

.855698   - 

•.921370   1.216290 

0 
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APPENDIX      C.  EXAMPLE  OUTPUT  DATA  FROM  PROGRAM 

USPOTF2 

DATA   READ   FROM   FILE   CODE    1 


XKXxwirwjr7rx«K«7r7r ^  x.  *.^ r.  m.  *.  *  ^  *  ^  ^  ^  m  f.  n  *  k  ^  j\  *  ^  *.  *.  *.  x.  w.  ^  jk  a  *,  a  m.  ^  ^ 71  ^  x  x.  *  ^  Jnnnonnnrmi 

AIRFOIL  TYPE    :    8.4X  THICKNESS  SYMMETRICAL  MISES  AIRFOIL 
AIRFOIL   COORDINATES  ARE   USER   INPUT   DATA   (IFLAG   =    1) 
NLOWER   =    25    ,   NUPPER   =   25 

1  25        25 

2  0.0      2.0 


1. 

000000 

0. 

994858 

0. 

980866 

0. 

958884 

0. 

929536 

0. 

893455 

0. 

851308 

0. 

803815 

0. 

751753 

0. 

695948 

0. 

637271 

0. 

576620 

0. 

514918 

0. 

453098 

0. 

392084 

0. 

332794 

0. 

276105 

0. 

222865 

0. 

173861 

0. 

129819 

0, 

091393 

0. 

059146 

0. 

033560 

0. 

015010 

0. 

003767 

0. 

000000 

0. 

003767 

0. 

015008 

0, 

033560 

0, 

059146 

0. 

091393 

0. 

129819 

0, 

173861 

0. 

222865 

0. 

276105 

0. 

332791 

0, 

392082 

0. 

453095 

0. 

514915 

0. 

576617 

0. 

637266 

0. 

695946 

0. 

751750 

0. 

803815 

0. 

851308 

0. 

893455 

0. 

929536 

0, 

958884 

0. 

980866 

0. 

994858 

1, 

000000 

0. 

000000 

-0. 

000782 

-0, 

002784 

-0, 

005721 

-0. 

009351 

-0, 

,013459 

-0. 

017837 

-0. 

022285 

-0, 

026618 

-0. 

030671 

-0. 

034289 

-0, 

,037341 

-0. 

039712 

-0. 

041314 

-0, 

,042083 

-0. 

041979 

-0. 

040979 

-0, 

,039096 

-0, 

036360 

-0. 

032820 

-0. 

,028555 

-0. 

023651 

-0, 

018220 

-0, 

,012379 

-0. 

.006259 

0. 

000000 

0, 

.006259 

0. 

012379 

0, 

018220 

0 

,023651 

0. 

028555 

0, 

032820 

0 

,036360 

0. 

039096 

0, 

,040979 

0, 

,041979 

0, 

,042083 

0, 

,041314 

0. 

,039712 

0. 

,037341 

0. 

,034269 

0. 

,030671 

0. 

,026618 

0. 

,022285 

0. 

.017837 

0 

,013459 

0 

,009351 

0. 

,005721 

0. 

,002784 

0, 

,000782 

0 

,000000 

1, 

,000000 

0, 

, 994858 

0. 

, 980866 

0 

, 958884 

0. 

,929536 

0 

,893455 

0 

,851308 

0. 

.803815 

0. 

,751753 

0 

.695948 

0 

,637271 

0 

,576620 

0 

.514918 

0 

.453098 

0 

,392084 

0 

.332794 

0 

,276105 

0 

.222865 

0 

.173861 

0, 

,129819 

0 

,091393 

0 

.059146 

0 

,033560 

0 

.015010 

0 

.003767 

0 

.000000 

0 

,003767 

0 

.015008 

0 

,033560 

0 

.059146 

0 

.091393 

0. 

.129819 

0 

,173861 

0 

,222865 

0 

,276105 

0 

.332791 

0 

.392082 

0 

,453095 

0 

.514915 

0 

,576617 

0 

,637266 

0 

.695946 

0 

,751750 

0 

,803815 

0 

.851308 

0 

,893455 

0 

.929536 

0 

, 958884 

0 

.980866 

0 

, 994858 

1 

.000000 

0 

.000000 

-0 

,000782 

-0 

.002784 

-0 

.005721 

-0 

.009351 

-0 

.013459 

-0 

.017837 

-0 

,022285 

-0 

.026618 

-0 

.030671 

-0 

.034289 

-0 

.037341 

-0 

.039712 

-0 

,041314 

-0 

.042083 

-0 

.041979 

-0 

.040979 

-0 

.039096 

-0 

.036360 

-0 

.032820 

-0 

.028555 

-0 

.023651 

-0 

.018220 

-0 

.012379 

-0 

.006259 

0 

.000000 

0 

.006259 

0 

.012379 

0. 

.018220 

0. 

.023651 

0 

.028555 

0 

.032820 

0 

.036360 

0 

.039096 

0 

.040979 

0 

.041979 

0 

.042083 

0 

.041314 

0 

.039712 

0 

.037341 

0 

.034289 

0 

.030671 

0 

.026618 

0 

.022285 

0 

.017837 

0 

.013459 

0 

.009351 

0. 

.005721 

0 

.002784 

0 

.000782 

0 

.000000 

0 

.000000 

0 

.000000- 

45 

.800003 

0 

.000000 

0 

.000000 

0 

.000000 

0.000000 

0 

0 

.000000 

0 

.000000 

0 

.000000 

0 

.000000 

0 

.000000 

0 

.000000 

0.000000 

0.000000 

0 

.650000 

0 

,025000 

0 

.001000 

0 

.000000 

4 

.855698 

-0 

.921370 

1.216290 

0 
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